DiggingCPC_4

Analyzing Multiple Plays of Digging Game Data

library(readr)
library(NatParksPalettes)
library(ggplot2)
library(Hmisc) #for computing correlation stuff

## 
## Attaching package: 'Hmisc'

## The following objects are masked from 'package:base':
## 
##     format.pval, units

library(knitr) 
library(tidyverse, warn.conflict=F)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ stringr   1.5.1
## ✔ forcats   1.0.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter()    masks stats::filter()
## ✖ dplyr::lag()       masks stats::lag()
## ✖ dplyr::src()       masks Hmisc::src()
## ✖ dplyr::summarize() masks Hmisc::summarize()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(corrplot)

## corrplot 0.95 loaded

library(dplyr)
library(MASS)

## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select

library(sandwich)
library(lmtest)

## Loading required package: zoo
## 
## Attaching package: 'zoo'
## 
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(broom)
library(dplyr)
library(gt)

## 
## Attaching package: 'gt'
## 
## The following object is masked from 'package:Hmisc':
## 
##     html

library(lme4)

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

library(dunn.test)

source("poster_theme.R") #for figure themes

## Warning: The `size` argument of `element_line()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

source("functions/process_initquiz.R")
source("functions/process_digging_long.R")
source("functions/process_params.R")
source("functions/softmax.R")

# RUN EVERYTIME YOU SWITHC PLAYS
#data_path <- "C:\\Users\\Ellen Martin\\OneDrive\\Desktop\\Rutledge Lab\\ApathyDepression\\updatedData\\data\\happyData\\"
data_path <- "~/Desktop/Rutledge Lab/ApathyDepression/updatedData/data/happyData/"

# play one for people with at least one play 
#one_long_file_z <- paste0(data_path,"z\\THP_uk_DiggingHappy2_z_n2214.csv") #zscored
one_long_file_z <- paste0(data_path,"z/THP_uk_DiggingHappy2_z_n2214.csv") #zscored

# parameters and model fits data (4 parameter model)
#happy_file_one_z <- paste0(data_path, "modelBased\\z\\happy_params4_play2_z_n2214.csv")
happy_file_one_z <- paste0(data_path, "modelBased/z/happy_params4_play2_z_n2214.csv")

# 6 parameter model
#happy_file_one_z_6 <- paste0(data_path, "modelBased\\z\\happy_params6_play2_z_n2214.csv")
happy_file_one_z_6 <- paste0(data_path, "modelBased/z/happy_params6_play2_z_n2214.csv")

# survey data
#survey_file <- paste0("C:\\Users\\Ellen Martin\\OneDrive\\Desktop\\Rutledge Lab\\ApathyDepression\\updatedData\\analyses\\uk_general_survey.csv")
survey_file <- paste0("~/Desktop/Rutledge Lab/ApathyDepression/updatedData/analyses/uk_general_survey.csv")

# initial quiz data
initquiz_file <- paste0(data_path, "THP_uk_initquiz.csv")

# Merging with survey data and backfilling


T_digging_long_z <- process_digging_long_z(one_long_file_z)

## Rows: 66416 Columns: 18
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): userKey
## dbl (17): TrialNumber, SafeProb, SafeValue, RiskyProb, RiskyValue, RiskySide...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

length(unique(T_digging_long_z$userKey)) #N=2200

## [1] 2200

T_digging_long_z <- process_params_z(T_digging_long_z,happy_file_one_z)

## Rows: 2214 Columns: 9
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): userKey
## dbl (8): b_evrpe_1, b_evrpe_2, b_evrpe_3, b_evrpe_4, r2, sse, aic, bic
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

length(unique(T_digging_long_z$userKey)) #2200

## [1] 2200

head(T_digging_long_z)

##        userKey  Variance TrialNumber SafeProb SafeValue RiskyProb RiskyValue
## 1 bemndakkrzxz 0.9999962           1      0.8       -20       0.1        -65
## 2 bemndakkrzxz 0.9999962           2      0.8       -20       0.3        -40
## 3 bemndakkrzxz 0.9999962           3      0.8       -20       0.4        -25
## 4 bemndakkrzxz 0.9999962           4      0.8       -20       0.7        -25
## 5 bemndakkrzxz 0.9999962           5      0.8        20       0.4         25
## 6 bemndakkrzxz 0.9999962           6      0.8        20       0.5         80
##   RiskySide Choice Outcome RevealDepth NextIsland LastIsland   zHappy HappyRT
## 1         1      2       0       0.286          2          1       NA      NA
## 2         1      2       0       1.100          2          1       NA      NA
## 3         1      2       0       0.429          2          1       NA      NA
## 4         1      2     -25       2.000          2          1 0.319487   1.669
## 5         2      2      25       2.000          1          1       NA      NA
## 6         2      2      80       2.000          1          1       NA      NA
##   PlayNo col15 col16 zHappyPred Trial ev_chosen_multi_z rpe_chosen_multi_z
## 1      2     0     0         NA  Loss      1.761076e-06         0.02451529
## 2      2     0     0         NA  Loss      1.761076e-06         0.02451529
## 3      2     0     0         NA  Loss      1.761076e-06         0.02451529
## 4      2     0     0  -0.511399  Loss      1.761076e-06         0.02451529
## 5      2     0     0         NA  Gain      1.761076e-06         0.02451529
## 6      2     0     0         NA  Gain      1.761076e-06         0.02451529
##   tau_multi_z const1_multi_z r2_multi_z sse_multi_z aic_multi_z bic_multi_z
## 1   0.2141355     -0.3948346  0.5768366    2.115817   -68.31738   -58.98856
## 2   0.2141355     -0.3948346  0.5768366    2.115817   -68.31738   -58.98856
## 3   0.2141355     -0.3948346  0.5768366    2.115817   -68.31738   -58.98856
## 4   0.2141355     -0.3948346  0.5768366    2.115817   -68.31738   -58.98856
## 5   0.2141355     -0.3948346  0.5768366    2.115817   -68.31738   -58.98856
## 6   0.2141355     -0.3948346  0.5768366    2.115817   -68.31738   -58.98856

happyData1 <- process_survey_data(T_digging_long_z,survey_file)

## Rows: 11142 Columns: 58
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): userKey
## dbl (57): phq8_total, phq8_q1, phq8_q2, phq8_q3, phq8_q4, phq8_q5, phq8_q6, ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

happyData1 <- happyData1 %>%
  arrange(userKey, PlayNo, TrialNumber)


# backfill z-scored
happyData1 <- happyData1 %>%
  mutate(zHappy_filled = zHappy) %>%
  fill(zHappy_filled, .direction = "up")

happyData1 <- happyData1 %>%
  # Backfill NaNs in zHappyRating and store in zHappy column
  mutate(zHappyPred_filled = zHappyPred) %>%
  fill(zHappyPred_filled, .direction = "up")

length(unique(happyData1$userKey))

## [1] 1922

# GAD and PHQ bins - may have to alter this based on medigan GAD in this play
happyData1 <- happyData1 %>%
  mutate(GAD_bin = cut(gad7_total, breaks = c(-Inf, 4, 9, 14, 21), labels = c("0-4", "5-9", "10-14", "15-21")),
      PHQ_bin = cut(phq8_total, breaks = c(-Inf, 4, 9, 14, 24), labels = c("0-4", "5-9", "10-14", "15-24")),
      GAD_binary = case_when(gad7_total >= 6 ~ "GAD 6+",
                             gad7_total < 6 ~ "GAD <6"),
      PHQ_binary = case_when(phq8_total >= 7 ~ "PHQ 7+",
                             phq8_total < 7 ~ "PHQ <7"))

# Anxiety and Depression Diagnoses Analyses
T_survey <- read_csv(survey_file)

## Rows: 11142 Columns: 58
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): userKey
## dbl (57): phq8_total, phq8_q1, phq8_q2, phq8_q3, phq8_q4, phq8_q5, phq8_q6, ...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

diag <- T_survey %>%
  dplyr::select(c("diagnosis_depression","diagnosis_anxiety","userKey"))

head(diag)

## # A tibble: 6 × 3
##   diagnosis_depression diagnosis_anxiety userKey     
##                  <dbl>             <dbl> <chr>       
## 1                  NaN               NaN abemndvlknzy
## 2                  NaN               NaN apwrxpqlwnmv
## 3                    0                 0 bemndakkrzxz
## 4                    0                 0 bemndbaadnzy
## 5                  NaN               NaN bemndbawdnzy
## 6                    0                 0 bemndbbrqnzy

happyData1_diag <- merge(happyData1, diag, by="userKey")

happyData1_diag <- happyData1_diag %>%
  filter(!is.na(diagnosis_anxiety) & !is.na(diagnosis_depression)) %>%
  mutate(diagnosis_anxiety = as.factor(diagnosis_anxiety),
         diagnosis_depression = as.factor(diagnosis_depression)) %>%
  mutate(diagnosis_anxiety = case_when(diagnosis_anxiety == 1 ~ "anx",
                                       diagnosis_anxiety == 0 ~ "no anx"),
         diagnosis_depression = case_when(diagnosis_depression == 1 ~ "dep",
                                          diagnosis_depression == 0 ~ "no dep"))


length(unique(happyData1_diag$userKey)) 

## [1] 1585

sample_size = length(unique(happyData1_diag$userKey))

# past
happyData1_past <- happyData1 %>%
  arrange(userKey, TrialNumber) %>%
  group_by(userKey) %>%
  mutate(
    # Create an indicator for valid happiness ratings
    happyind = !is.na(zHappy),
    
    # Identify if the current trial is on the same island as the previous valid happiness rating
    prev_happyind_trial = lag(Trial, order_by = TrialNumber),  # Previous trial type (Gain/Loss)
    
    # Set PastIsland based on the previous valid happiness rating's trial type
    PastIsland = case_when(
      lag(happyind, order_by = TrialNumber) & prev_happyind_trial == "Gain" ~ 2,
      lag(happyind, order_by = TrialNumber) & prev_happyind_trial == "Loss" ~ 1,
      TRUE ~ NA_real_  # NA for the first island or no previous valid row
    )
  ) %>%
  ungroup()

happyData1_past <- happyData1_past %>%
  fill(PastIsland, .direction = "up")

sample_size = length(unique(happyData1$userKey))

saveRDS(happyData1, "~/Desktop/Rutledge Lab/ApathyDepression/updatedData/data/happyData/happyData1.rds")
saveRDS(happyData1_diag, "~/Desktop/Rutledge Lab/ApathyDepression/updatedData/data/happyData/happyData1_diag.rds")
saveRDS(happyData1_past, "~/Desktop/Rutledge Lab/ApathyDepression/updatedData/data/happyData/happyData1_past.rds")

Reading in 6 Parameter Happiness Model

Read this in too and then return to main script DiggingCPC_5_play1.

params6 <- read.csv(happy_file_one_z_6)
params_multi <- params6 %>%
    rename(
      ev_safe = b_evrpe_1,
      rpe_safe = b_evrpe_2,
      ev_risky = b_evrpe_3,
      rpe_risky = b_evrpe_4,
      tau_split = b_evrpe_5,
      const_split = b_evrpe_6,
      r2_split = r2,
      sse_split = sse,
      aic_split = aic,
      bic_split =bic
    )

happyData1_6 <- merge(happyData1, params_multi, by="userKey")
mean(params_multi$r2_split)

## [1] 0.9375163

mean(params_multi$aic_split)

## [1] -196.5986

mean(params_multi$bic_split)

## [1] -182.6054

# mean tau and baseline happiness

mean(params_multi$tau_split) 

## [1] 0.5539281

mean(params_multi$const_split)

## [1] -0.0003904446

cor.test(params_multi$ev_risky,params_multi$ev_safe,method="spearman")

## Warning in cor.test.default(params_multi$ev_risky, params_multi$ev_safe, :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  params_multi$ev_risky and params_multi$ev_safe
## S = 2534661065, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.4013234

cor.test(params_multi$rpe_risky,params_multi$rpe_safe,method="spearman")

## Warning in cor.test.default(params_multi$rpe_risky, params_multi$rpe_safe, :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  params_multi$rpe_risky and params_multi$rpe_safe
## S = 1800574074, p-value = 0.8314
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## 0.004527008

cor.test(params_multi$rpe_risky,params_multi$ev_risky,method="spearman")

## Warning in cor.test.default(params_multi$rpe_risky, params_multi$ev_risky, :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  params_multi$rpe_risky and params_multi$ev_risky
## S = 1800296266, p-value = 0.8258
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## 0.004680598

cor.test(params_multi$rpe_safe,params_multi$ev_safe,method="spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  params_multi$rpe_safe and params_multi$ev_safe
## S = 2198357878, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.2153934

% of gambling

risky_data <- happyData1 %>%
    group_by(userKey, Trial) %>%
    summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop')


summary_data <- risky_data %>%
    group_by(Trial) %>%
    summarise(mean_percent_risky = mean(percent_risky, na.rm = TRUE),
              se_percent_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
              .groups = 'drop')
  
  ggplot(summary_data, aes(x = Trial, y = mean_percent_risky, color = Trial)) +
    geom_point(position = position_dodge(0.8), size = 3) +
    geom_errorbar(aes(ymin = mean_percent_risky - se_percent_risky, ymax = mean_percent_risky + se_percent_risky),
                  width = 0.2) +
    labs(
      title = paste0("% Risky Choices (Play 2) (N= ", sample_size, ")"),
      x = "Trial Type",
      y = "% Risky Choices"
    ) +
    theme_minimal() +
    ylim(40, 60) +  
    poster_theme + 
    scale_color_manual(values = c("Gain" = "gold", "Loss" = "coral"))

- increased gambling in gain vs loss trials

z-scored happiness

#no correlation between PHQ or GAD with z-scored happiness in gain or loss domain
summary_zHappy_1 <- happyData1 %>%
    group_by(userKey) %>%
    summarise(
      mean_zHappy_filled = mean(zHappy_filled, na.rm=TRUE),
      sd_zHappy_filled = sd(zHappy_filled, na.rm = TRUE),
      rsd_Happy_filled = relSD_tc(zHappy, -2, 2),
      GAD_score = first(gad7_total),
      PHQ_score = first(phq8_total),
      GAD_bin = first(GAD_bin),
      PHQ_bin = first(PHQ_bin),
      gamble_decisions = sum(Choice == 2, na.rm = TRUE),
      total_decisions = n(),
      percent_gamble = (gamble_decisions / total_decisions) * 100,
      const = first(const1_multi_z),
      tau = first(tau_multi_z),
      rpe = first(rpe_chosen_multi_z),
      ev = first(ev_chosen_multi_z)
    ) 

cor.test(summary_zHappy_1$sd_zHappy_filled,summary_zHappy_1$GAD_score,method="spearman") # near sig correlation between GAD and sd zHappy

## Warning in cor.test.default(summary_zHappy_1$sd_zHappy_filled,
## summary_zHappy_1$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  summary_zHappy_1$sd_zHappy_filled and summary_zHappy_1$GAD_score
## S = 1168401382, p-value = 0.54
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.01404795

cor.test(summary_zHappy_1$sd_zHappy_filled,summary_zHappy_1$PHQ_score,method="spearman") # not sig with PHQ

## Warning in cor.test.default(summary_zHappy_1$sd_zHappy_filled,
## summary_zHappy_1$PHQ_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  summary_zHappy_1$sd_zHappy_filled and summary_zHappy_1$PHQ_score
## S = 1213557765, p-value = 0.2631
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## -0.0255378

cor.test(summary_zHappy_1$rsd_Happy_filled,summary_zHappy_1$GAD_score,method="spearman") # no correlation between GAD and relSD

## Warning in cor.test.default(summary_zHappy_1$rsd_Happy_filled,
## summary_zHappy_1$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  summary_zHappy_1$rsd_Happy_filled and summary_zHappy_1$GAD_score
## S = 1121343489, p-value = 0.2425
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.02679329

cor.test(summary_zHappy_1$rsd_Happy_filled,summary_zHappy_1$PHQ_score,method="spearman") # no correlation between PHQ and relSD

## Warning in cor.test.default(summary_zHappy_1$rsd_Happy_filled,
## summary_zHappy_1$PHQ_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  summary_zHappy_1$rsd_Happy_filled and summary_zHappy_1$PHQ_score
## S = 1164837936, p-value = 0.4934
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.01563373

#GAD and sd Happiness (z-scored)
summary_bins <- summary_zHappy_1 %>%
    group_by(GAD_bin) %>%
    summarise(
      mean_sd_Happy_filled = mean(mean_zHappy_filled, na.rm = TRUE),
      sem_sd_Happy_filled = sd(mean_zHappy_filled, na.rm = TRUE) / sqrt(n())
    ) %>% na.omit()

ggplot(summary_bins, aes(x = GAD_bin, y = mean_sd_Happy_filled)) +
    geom_point(size = 3) +
    geom_errorbar(aes(ymin = mean_sd_Happy_filled - sem_sd_Happy_filled, ymax = mean_sd_Happy_filled + sem_sd_Happy_filled),
                  width = 0.2) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    labs(
      title = "Happiness by Anxiety (Play 2) N=1922",
      x = "Anxiety (GAD score)",
      y = "z-scored Happiness"
    ) +
    theme_minimal() +
    ylim(-0.02, 0.02) +
    poster_theme 

# by Domain
summary_zHappy_1 <- happyData1 %>%
    group_by(userKey,Trial) %>%
    summarise(
      mean_zHappy_filled = mean(zHappy_filled, na.rm=TRUE),
      sd_zHappy_filled = sd(zHappy_filled, na.rm = TRUE),
      rsd_Happy_filled = relSD_tc(zHappy, -2, 2),
      GAD_score = first(gad7_total),
      PHQ_score = first(phq8_total),
      GAD_bin = first(GAD_bin),
      PHQ_bin = first(PHQ_bin),
      gamble_decisions = sum(Choice == 2, na.rm = TRUE),
      total_decisions = n(),
      percent_gamble = (gamble_decisions / total_decisions) * 100,
      const = first(const1_multi_z),
      tau = first(tau_multi_z),
      rpe = first(rpe_chosen_multi_z),
      ev = first(ev_chosen_multi_z)
    ) 

## `summarise()` has grouped output by 'userKey'. You can override using the
## `.groups` argument.

summary_bins <- summary_zHappy_1 %>%
    group_by(GAD_bin,Trial) %>%
    summarise(
      mean_z_Happy_filled = mean(mean_zHappy_filled, na.rm = TRUE),
      sem_sd_Happy_filled = sd(mean_zHappy_filled, na.rm = TRUE) / sqrt(n())
    ) %>% na.omit()

## `summarise()` has grouped output by 'GAD_bin'. You can override using the
## `.groups` argument.

ggplot(summary_bins, aes(x = GAD_bin, y = mean_z_Happy_filled, color=Trial)) +
    geom_point(size = 3) +
    geom_errorbar(aes(ymin = mean_z_Happy_filled - sem_sd_Happy_filled, ymax = mean_z_Happy_filled + sem_sd_Happy_filled),
                  width = 0.2) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    labs(
      title = "Happiness by Anxiety",
      x = "Anxiety (GAD score)",
      y = "z-scored Happiness"
    ) +
  scale_color_manual(values = c("Gain" = "gold", "Loss" = "coral")) + 
    theme_minimal() +
    ylim(-0.2, 0.2)+
    poster_theme 

Happiness Model Parameters

summary_params <- happyData1 %>%
    dplyr::select(ev_chosen_multi_z, rpe_chosen_multi_z) %>%
    pivot_longer(cols = c(ev_chosen_multi_z, rpe_chosen_multi_z), 
                 names_to = "parameter", 
                 values_to = "value") %>%
    group_by(parameter) %>%
    summarise(
      mean_value = mean(value, na.rm = TRUE),
      se_value = sd(value, na.rm = TRUE) / sqrt(n())
    ) %>%
    na.omit()
  
  # Generate the plot
  ggplot(summary_params, aes(x = parameter, y = mean_value)) +
    geom_bar(stat = "identity", position = position_dodge(), width = 0.7, fill = "lightblue", color = "black") +
    geom_errorbar(aes(ymin = mean_value - se_value, ymax = mean_value + se_value),
                  position = position_dodge(0.7), width = 0.2) +
    geom_hline(yintercept = 0, linetype = "solid", color = "black", size = 0.5) +
    labs(
      title = "Play 2 Parameters (N=1922)",
      x = "Parameter",
      y = "Parameter Estimate"
    ) +
    theme_minimal() +
   poster_theme  + 
    scale_x_discrete(labels = c(
      "ev_chosen_multi_z" = "EV", 
      "rpe_chosen_multi_z" = "RPE"))

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# split by GAD
summary_params <- happyData1 %>%
  dplyr::select(GAD_binary, ev_chosen_multi_z, rpe_chosen_multi_z) %>%
  group_by(GAD_binary) %>% 
    pivot_longer(cols = c(ev_chosen_multi_z, rpe_chosen_multi_z), 
                 names_to = "parameter", 
                 values_to = "value") %>%
    group_by(parameter, GAD_binary) %>%
    summarise(
      mean_value = mean(value, na.rm = TRUE),
      se_value = sd(value, na.rm = TRUE) / sqrt(n())
    ) %>%
    na.omit() 

## `summarise()` has grouped output by 'parameter'. You can override using the
## `.groups` argument.

  # Generate the plot
  ggplot(summary_params, aes(x = parameter, y = mean_value, fill = GAD_binary)) +
    geom_bar(stat = "identity", position = position_dodge(), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = mean_value - se_value, ymax = mean_value + se_value),
                  position = position_dodge(0.7), width = 0.2) +
    labs(
      x = "Parameter",
      y = "Parameter Estimate",
    ) +
    theme_minimal() +
    poster_theme + 
    scale_fill_manual(values = c("GAD 6+" = "#00BD9D", "GAD <6" = "#8BD7D2")) +
    scale_x_discrete(labels = c("ev_chosen_multi_z" = "EV", "rpe_chosen_multi_z" = "RPE"))

  # split by PHQ
  summary_params <- happyData1 %>%
  dplyr::select(PHQ_binary, ev_chosen_multi_z, rpe_chosen_multi_z) %>%
  group_by(PHQ_binary) %>% 
    pivot_longer(cols = c(ev_chosen_multi_z, rpe_chosen_multi_z), 
                 names_to = "parameter", 
                 values_to = "value") %>%
    group_by(parameter, PHQ_binary) %>%
    summarise(
      mean_value = mean(value, na.rm = TRUE),
      se_value = sd(value, na.rm = TRUE) / sqrt(n())
    ) %>%
    na.omit() 

## `summarise()` has grouped output by 'parameter'. You can override using the
## `.groups` argument.

  # Generate the plot
  ggplot(summary_params, aes(x = parameter, y = mean_value, fill = PHQ_binary)) +
    geom_bar(stat = "identity", position = position_dodge(), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = mean_value - se_value, ymax = mean_value + se_value),
                  position = position_dodge(0.7), width = 0.2) +
    labs(
      x = "Parameter",
      y = "Parameter Estimate",
    ) +
    theme_minimal() +
    poster_theme + 
    scale_fill_manual(values = c("PHQ 7+" = "red", "PHQ <7" = "pink")) +
    scale_x_discrete(labels = c("ev_chosen_multi_z" = "EV", "rpe_chosen_multi_z" = "RPE"))

## correlations between parameters and symptoms

cor.test(summary_zHappy_1$GAD_score,summary_zHappy_1$rpe,method="spearman") #significant negative correlation between anxiety and rpe

## Warning in cor.test.default(summary_zHappy_1$GAD_score, summary_zHappy_1$rpe, :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  summary_zHappy_1$GAD_score and summary_zHappy_1$rpe
## S = 9187684263, p-value = 0.8406
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## 0.003258788

cor.test(summary_zHappy_1$GAD_score,summary_zHappy_1$ev,method="spearman") #n.s. corr between anxiety and ev

## Warning in cor.test.default(summary_zHappy_1$GAD_score, summary_zHappy_1$ev, :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  summary_zHappy_1$GAD_score and summary_zHappy_1$ev
## S = 8757436081, p-value = 0.002048
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.04993498

cor.test(summary_zHappy_1$PHQ_score,summary_zHappy_1$rpe,method="spearman") #more significant negative correlation betwen depression and rpe

## Warning in cor.test.default(summary_zHappy_1$PHQ_score, summary_zHappy_1$rpe, :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  summary_zHappy_1$PHQ_score and summary_zHappy_1$rpe
## S = 9290631743, p-value = 0.249
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.01859924

cor.test(summary_zHappy_1$PHQ_score,summary_zHappy_1$ev,method="spearman") #n.s. corr between depression and ev

## Warning in cor.test.default(summary_zHappy_1$PHQ_score, summary_zHappy_1$ev, :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  summary_zHappy_1$PHQ_score and summary_zHappy_1$ev
## S = 9238539054, p-value = 0.1352
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.02410197

• there is a more significant correlation between depression and rpe than anxiety and rpe in play 1

Winning vs Losing Happiness (Taking a Risk in Gain and Loss Domain)

result <- winvloss_z_happy(happyData1, "Gain", "Risky")
result$summary_long

## # A tibble: 2 × 3
##   OutcomeType mean_happiness se_happiness
##   <chr>                <dbl>        <dbl>
## 1 Loss                -0.322       0.0333
## 2 Win                  0.511       0.0308

result$plot

result$wilcox_test 

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  merged_data$Win and merged_data$Loss
## V = 128082, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

# model underestimates unhappiness when losing a gamble in the gain domain

result <- winvloss_z_happy_GAD(happyData1, "Gain", "Risky")

## Warning in cor.test.default(filtered_data$zHappy[filtered_data$OutcomeType == :
## Cannot compute exact p-value with ties

## Warning in cor.test.default(filtered_data$zHappy[filtered_data$OutcomeType == :
## Cannot compute exact p-value with ties

result$summary_long

## # A tibble: 8 × 4
##   OutcomeType GAD_bin mean_happiness se_happiness
##   <chr>       <fct>            <dbl>        <dbl>
## 1 Loss        0-4             -0.352       0.0505
## 2 Loss        5-9             -0.343       0.0611
## 3 Loss        10-14           -0.181       0.0811
## 4 Loss        15-21           -0.330       0.0996
## 5 Win         0-4              0.507       0.0470
## 6 Win         5-9              0.499       0.0586
## 7 Win         10-14            0.550       0.0767
## 8 Win         15-21            0.510       0.0853

result$plot

result$wilcox_test_gad6plus

## NULL

result$wilcox_test_gad_less6 

## NULL

result$mw_test_win

## NULL

result$mw_test_loss #GAD groups significantly different in happiness when losing a risky choice in loss domain

## NULL

#corelations
result$spearman_corr_win

## 
##  Spearman's rank correlation rho
## 
## data:  filtered_data$zHappy[filtered_data$OutcomeType == "Win"] and filtered_data$gad7_total[filtered_data$OutcomeType == "Win"]
## S = 335559632, p-value = 0.4388
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.02172739

result$spearman_corr_loss

## 
##  Spearman's rank correlation rho
## 
## data:  filtered_data$zHappy[filtered_data$OutcomeType == "Loss"] and filtered_data$gad7_total[filtered_data$OutcomeType == "Loss"]
## S = 394954110, p-value = 0.1761
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.03684359

# model is off for loss domain

#result <- winvloss_z_happy_resid_GAD(happyData1, "Gain", "Risky")
#result$summary_long
#result$plot
#result$wilcox_test_gad6plus
#result$wilcox_test_gad_less6 
#result$mw_test_win
#result$mw_test_loss

# model is overall overestimating happiness (since residual is negative, it means predicted > actual)

result <- winvloss_z_happy(happyData1, "Loss", "Risky")
result$summary_long

## # A tibble: 2 × 3
##   OutcomeType mean_happiness se_happiness
##   <chr>                <dbl>        <dbl>
## 1 Loss                -0.490       0.0305
## 2 Win                  0.233       0.0276

result$plot

result$wilcox_test 

## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  merged_data$Win and merged_data$Loss
## V = 167670, p-value < 2.2e-16
## alternative hypothesis: true location shift is not equal to 0

# model underestimates unhappiness when losing a gamble in the gain domain

result <- winvloss_z_happy_GAD(happyData1, "Loss", "Risky")

## Warning in cor.test.default(filtered_data$zHappy[filtered_data$OutcomeType == :
## Cannot compute exact p-value with ties

## Warning in cor.test.default(filtered_data$zHappy[filtered_data$OutcomeType == :
## Cannot compute exact p-value with ties

result$summary_long

## # A tibble: 8 × 4
##   OutcomeType GAD_bin mean_happiness se_happiness
##   <chr>       <fct>            <dbl>        <dbl>
## 1 Loss        0-4             -0.477       0.0486
## 2 Loss        5-9             -0.540       0.0537
## 3 Loss        10-14           -0.429       0.0819
## 4 Loss        15-21           -0.490       0.0799
## 5 Win         0-4              0.206       0.0452
## 6 Win         5-9              0.291       0.0507
## 7 Win         10-14            0.238       0.0633
## 8 Win         15-21            0.187       0.0687

result$plot

result$wilcox_test_gad6plus

## NULL

result$wilcox_test_gad_less6 

## NULL

result$mw_test_win

## NULL

result$mw_test_loss #GAD groups significantly different in happiness when winning a risky choice in loss domain - high anx subjects are less happy after winning a risky choice in loss domain

## NULL

#corelations
result$spearman_corr_win

## 
##  Spearman's rank correlation rho
## 
## data:  filtered_data$zHappy[filtered_data$OutcomeType == "Win"] and filtered_data$gad7_total[filtered_data$OutcomeType == "Win"]
## S = 1300621466, p-value = 0.07984
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.03959509

result$spearman_corr_loss

## 
##  Spearman's rank correlation rho
## 
## data:  filtered_data$zHappy[filtered_data$OutcomeType == "Loss"] and filtered_data$gad7_total[filtered_data$OutcomeType == "Loss"]
## S = 254913169, p-value = 0.5832
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.01622529

# model is off for loss domain

#result <- winvloss_z_happy_resid_GAD(happyData1, "Loss", "Risky")
#result$summary_long
#result$plot
#result$wilcox_test_gad6plus
#result$wilcox_test_gad_less6 
#result$mw_test_win
#result$mw_test_loss

# in loss domain, model is overestimating happiness, particularly for low anxiety participants when they win the gamble

Future Information

• Happiness based just on the next island
• Split by GAD, then GAD bins
• Split by PHQ, then PHQ bis

# happiness based on next island (positive or negative)


positive_happy <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2) %>%
  group_by(userKey) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE)) %>%
  mutate(Future = "Positive")
  
  # Calculate mean happiness and predicted happiness for NextIsland == 1 within subjects
negative_happy <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1) %>%
  group_by(userKey) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), .groups = 'drop') %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE)) %>%
  mutate(Future = "Negative")
  
  # Combine results
combined_happy <- bind_rows(positive_happy, negative_happy)
  
  # Plotting
  ggplot(combined_happy, aes(x = Future, y = overall_mean_happiness, fill = Future)) +
    geom_bar(stat = "identity", position = position_dodge(), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = overall_mean_happiness - sem_happiness, ymax = overall_mean_happiness + sem_happiness),
                  position = position_dodge(0.7), width = 0.2) +
    geom_point(aes(y = overall_mean_pred_happiness), color = "lightblue", shape = 8, size = 3) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0

    scale_fill_manual(values = c("Positive" = "#FFD700", "Negative" = "#FF6347")) +
    labs(
      title = "z-scored Happiness based on Future",
      x = "Future",
      y = "z-scored Happiness"
    ) +
    theme_minimal() +
    ylim(-0.1, 0.1) +
    poster_theme  

  # happier when future is negative
  
# residuals
  positive_happy_resid <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2) %>%
  group_by(userKey) %>%
  summarise(mean_residual = mean(zHappy - zHappyPred, na.rm = TRUE), .groups = 'drop') %>%
    summarise(overall_mean_residual = mean(mean_residual, na.rm = TRUE),
              sem_residual = sd(mean_residual, na.rm = TRUE) / sqrt(n())) %>%
  mutate(Future = "Positive")
  
  # Calculate mean happiness and predicted happiness for NextIsland == 1 within subjects
negative_happy_resid <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1) %>%
  group_by(userKey) %>%
  summarise(mean_residual = mean(zHappy - zHappyPred, na.rm = TRUE), .groups = 'drop') %>%
    summarise(overall_mean_residual = mean(mean_residual, na.rm = TRUE),
              sem_residual = sd(mean_residual, na.rm = TRUE) / sqrt(n())) %>%
  mutate(Future = "Negative")
  
  # Combine results
combined_happy <- bind_rows(positive_happy_resid, negative_happy_resid)
  
  # Plotting
  ggplot(combined_happy, aes(x = Future, y = overall_mean_residual, color = Future)) +
    geom_point(position = position_dodge(0.7), size = 1) +  # Points instead of bars
    geom_errorbar(aes(ymin = overall_mean_residual - sem_residual, ymax = overall_mean_residual + sem_residual),
                  position = position_dodge(0.7), width = 0.2) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0

    scale_color_manual(values = c("Positive" = "#FFD700", "Negative" = "#FF6347")) +
    labs(
      title = "Residual Happiness based on Future",
      x = "Future",
      y = "Resid. Happiness"
    ) +
    theme_minimal() +
    poster_theme  

## GAD
  
  positive_happy <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2) %>%
  group_by(userKey, GAD_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop') %>%
  mutate(Future = "Positive")

# Calculate mean happiness and predicted happiness for negative future (NextIsland == 1, LastIsland == 1)
negative_happy <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1) %>%
  group_by(userKey, GAD_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop') %>%
  mutate(Future = "Negative")

# Combine results
combined_happy <- bind_rows(positive_happy, negative_happy)
combined_happy <- na.omit(combined_happy)

# Plotting
ggplot(combined_happy, aes(x = Future, y = overall_mean_happiness, fill = GAD_binary, group = GAD_binary)) +
  geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
  geom_errorbar(aes(ymin = overall_mean_happiness - sem_happiness, ymax = overall_mean_happiness + sem_happiness),
                position = position_dodge(0.7), width = 0.2) +
  geom_point(aes(y = overall_mean_pred_happiness), color = "lightblue", shape = 8, size = 3, 
             position = position_dodge(0.7)) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
  scale_fill_manual(values = c("GAD 6+" = "#00BD9D", "GAD <6" = "#8BD7D2")) +
  labs(
    title = "z-scored Happiness based on Future",
    x = "Future",
    y = "z-scored Happiness"
  ) +
  theme_minimal() +
  poster_theme 

## GAD_bin

 positive_happy <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2) %>%
  group_by(userKey,GAD_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop') %>%
  mutate(Future = "Positive")

# Calculate mean happiness and predicted happiness for negative future (NextIsland == 1, LastIsland == 1)
negative_happy <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1) %>%
  group_by(userKey, GAD_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop') %>%
  mutate(Future = "Negative")

# Combine results
combined_happy <- bind_rows(positive_happy, negative_happy)
combined_happy <- na.omit(combined_happy)

# Plotting
ggplot(combined_happy, aes(x = Future, y = overall_mean_happiness, fill = GAD_bin, group = GAD_bin)) +
  geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
  geom_errorbar(aes(ymin = overall_mean_happiness - sem_happiness, ymax = overall_mean_happiness + sem_happiness),
                position = position_dodge(0.7), width = 0.2) +
  geom_point(aes(y = overall_mean_pred_happiness), color = "lightblue", shape = 8, size = 3, 
             position = position_dodge(0.7)) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    scale_fill_manual(values = c("0-4" = "#8BD7D2", "5-9" = "#6FACA8", "10-14" = "#53817E", "15-21" = "#385654")) +
  labs(
    title = "z-scored Happiness based on Future",
    x = "Future",
    y = "z-scored Happiness"
  ) +
  theme_minimal() +
  poster_theme

## PHQ BINARY
positive_happy <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2) %>%
  group_by(userKey, PHQ_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop') %>%
  mutate(Future = "Positive")

# Calculate mean happiness and predicted happiness for negative future (NextIsland == 1, LastIsland == 1)
negative_happy <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1) %>%
  group_by(userKey, PHQ_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop') %>%
  mutate(Future = "Negative")

# Combine results
combined_happy <- bind_rows(positive_happy, negative_happy)
combined_happy <- na.omit(combined_happy)

# Plotting
ggplot(combined_happy, aes(x = Future, y = overall_mean_happiness, fill = PHQ_binary, group = PHQ_binary)) +
  geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
  geom_errorbar(aes(ymin = overall_mean_happiness - sem_happiness, ymax = overall_mean_happiness + sem_happiness),
                position = position_dodge(0.7), width = 0.2) +
  geom_point(aes(y = overall_mean_pred_happiness), color = "lightblue", shape = 8, size = 3, 
             position = position_dodge(0.7)) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
  scale_fill_manual(values = c("PHQ 7+" = "red", "PHQ <7" = "pink")) +
  labs(
    title = "z-scored Happiness based on Future",
    x = "Future",
    y = "z-scored Happiness"
  ) +
  theme_minimal() +
  poster_theme 

# PHQ Bin
## GAD_bin

 positive_happy <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2) %>%
  group_by(userKey,PHQ_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop') %>%
  mutate(Future = "Positive")

# Calculate mean happiness and predicted happiness for negative future (NextIsland == 1, LastIsland == 1)
negative_happy <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1) %>%
  group_by(userKey, PHQ_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop') %>%
  mutate(Future = "Negative")

# Combine results
combined_happy <- bind_rows(positive_happy, negative_happy)
combined_happy <- na.omit(combined_happy)

# Plotting
ggplot(combined_happy, aes(x = Future, y = overall_mean_happiness, fill = PHQ_bin, group = PHQ_bin)) +
  geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
  geom_errorbar(aes(ymin = overall_mean_happiness - sem_happiness, ymax = overall_mean_happiness + sem_happiness),
                position = position_dodge(0.7), width = 0.2) +
  geom_point(aes(y = overall_mean_pred_happiness), color = "lightblue", shape = 8, size = 3, 
             position = position_dodge(0.7)) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
      scale_fill_manual(values = c("0-4" = "#FFC0CB", "5-9" = "#CC9AA2", "10-14" = "#99737A", "15-24" = "#664D51")) +

  labs(
    title = "z-scored Happiness based on Future",
    x = "Future",
    y = "z-scored Happiness"
  ) +
  theme_minimal() +
  poster_theme

Risky Choices Future Information

positive_risky <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2) %>%
  group_by(userKey) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  mutate(Future = "Positive")
  
  # Calculate mean happiness and predicted happiness for NextIsland == 1 within subjects
negative_risky <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1) %>%
  group_by(userKey) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  mutate(Future = "Negative")
  
  # Combine results
combined_risky <- bind_rows(positive_risky, negative_risky)
summary_combined <- combined_risky %>%
  group_by(Future) %>%
  summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()))

  # Plotting
  ggplot(summary_combined, aes(x = Future, y = mean_risky, color = Future)) +
    geom_point(position = position_dodge(0.7), width = 0.7) +
    geom_errorbar(aes(ymin = mean_risky - sem_risky, ymax = mean_risky + sem_risky),
                  position = position_dodge(0.7), width = 0.2) + 
    geom_hline(yintercept = 50, linetype = "dotted", color = "black") +  # Add a dotted line at y=0

    scale_color_manual(values = c("Positive" = "#FFD700", "Negative" = "#FF6347")) +
    labs(
      title = "% Risky Choices based on Future (Play 2) N=1922",
      x = "Future",
      y = "% Risky Choices"
    ) +
    ylim(40,70) + 
    theme_minimal() +
    poster_theme 

## Warning in geom_point(position = position_dodge(0.7), width = 0.7): Ignoring
## unknown parameters: `width`

  # people choose to gamble more when the future is positive 
  
  #GAD
  positive_risky <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2) %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
    group_by(GAD_binary) %>%
    summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Positive")
  
  # non-paramateric t-test
  positive_risky_sig <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2) %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
            GAD_score = first(gad7_total),
              .groups = 'drop')
  
  wilcox.test(positive_risky_sig$percent_risky~positive_risky_sig$GAD_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  positive_risky_sig$percent_risky by positive_risky_sig$GAD_binary
## W = 455269, p-value = 0.8908
## alternative hypothesis: true location shift is not equal to 0

  cor.test(positive_risky_sig$percent_risky,positive_risky_sig$GAD_score,method="spearman")

## Warning in cor.test.default(positive_risky_sig$percent_risky,
## positive_risky_sig$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  positive_risky_sig$percent_risky and positive_risky_sig$GAD_score
## S = 1142249563, p-value = 0.706
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## 0.008649042

# no significant difference between GAD groups in risk taking when future is positive
  
  # Calculate mean happiness and predicted happiness for NextIsland == 1 within subjects
negative_risky <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1) %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  group_by(GAD_binary) %>%
  summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Negative")

negative_risky_sig <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1) %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
            GAD_score = first(gad7_total),
              .groups = 'drop')
  
  wilcox.test(negative_risky_sig$percent_risky~negative_risky_sig$GAD_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  negative_risky_sig$percent_risky by negative_risky_sig$GAD_binary
## W = 452774, p-value = 0.9432
## alternative hypothesis: true location shift is not equal to 0

    cor.test(negative_risky_sig$percent_risky,negative_risky_sig$GAD_score,method="spearman")

## Warning in cor.test.default(negative_risky_sig$percent_risky,
## negative_risky_sig$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  negative_risky_sig$percent_risky and negative_risky_sig$GAD_score
## S = 1.129e+09, p-value = 0.3797
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.020139

  # Combine results
combined_risky <- bind_rows(positive_risky, negative_risky)
combined_risky <- na.omit(combined_risky)

  # Plotting
  ggplot(combined_risky, aes(x = Future, y = mean_risky,  color = GAD_binary, group = GAD_binary)) +
    geom_point(position = position_dodge(0.7), width = 0.7) +
    geom_errorbar(aes(ymin = mean_risky - sem_risky, ymax = mean_risky + sem_risky),
                  position = position_dodge(0.7), width = 0.2) + 
    geom_hline(yintercept = 50, linetype = "dotted", color = "black") +  # Add a dotted line at y=0

    scale_color_manual(values = c("GAD 6+" = "#00BD9D", "GAD <6" = "#8BD7D2")) +
    labs(
      title = "% Risky Choices based on Future (Play 2) N=1922",
      x = "Future",
      y = "% Risky Choices"
    ) +
    ylim(40,70) + 
    theme_minimal() +
    poster_theme 

## Warning in geom_point(position = position_dodge(0.7), width = 0.7): Ignoring
## unknown parameters: `width`

    # GAD Bin
  
 positive_risky <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & !is.na(GAD_bin)) %>%
  group_by(userKey,GAD_bin) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
    group_by(GAD_bin) %>%
    summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Positive")
    
    
  
  negative_risky <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & !is.na(GAD_bin)) %>%
  group_by(userKey,GAD_bin) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  group_by(GAD_bin) %>%
  summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Negative")


  # Combine results
combined_risky <- bind_rows(positive_risky, negative_risky)
combined_risky <- na.omit(combined_risky)

  # Plotting
  ggplot(combined_risky, aes(x = Future, y = mean_risky,  color = GAD_bin, group = GAD_bin)) +
    geom_point(position = position_dodge(0.7), width = 0.7) +
    geom_errorbar(aes(ymin = mean_risky - sem_risky, ymax = mean_risky + sem_risky),
                  position = position_dodge(0.7), width = 0.2) + 
    geom_hline(yintercept = 50, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    scale_color_manual(values = c("0-4" = "#8BD7D2", "5-9" = "#6FACA8", "10-14" = "#53817E", "15-21" = "#385654")) +
    labs(
      title = "% Risky Choices based on Future (Play 4) N=1922",
      x = "Future",
      y = "% Risky Choices"
    ) +
        ylim(40,70) +
    theme_minimal() +
    poster_theme 

## Warning in geom_point(position = position_dodge(0.7), width = 0.7): Ignoring
## unknown parameters: `width`

  # PHQ binary
  
  positive_risky <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2) %>%
  group_by(userKey,PHQ_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
    group_by(PHQ_binary) %>%
    summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Positive")
  
  # non-paramateric t-test
  positive_risky_sig <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2) %>%
  group_by(userKey,PHQ_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
            PHQ_score = first(phq8_total),
              .groups = 'drop')
  
  wilcox.test(positive_risky_sig$percent_risky~positive_risky_sig$PHQ_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  positive_risky_sig$percent_risky by positive_risky_sig$PHQ_binary
## W = 465320, p-value = 0.7442
## alternative hypothesis: true location shift is not equal to 0

  cor.test(positive_risky_sig$percent_risky,positive_risky_sig$PHQ_score,method="spearman")

## Warning in cor.test.default(positive_risky_sig$percent_risky,
## positive_risky_sig$PHQ_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  positive_risky_sig$percent_risky and positive_risky_sig$PHQ_score
## S = 1169753384, p-value = 0.615
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.01147985

  # Calculate mean happiness and predicted happiness for NextIsland == 1 within subjects
negative_risky <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1) %>%
  group_by(userKey,PHQ_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  group_by(PHQ_binary) %>%
  summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Negative")

negative_risky_sig <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1) %>%
  group_by(userKey,PHQ_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
            PHQ_score = first(phq8_total),
              .groups = 'drop')
  
  wilcox.test(negative_risky_sig$percent_risky~negative_risky_sig$PHQ_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  negative_risky_sig$percent_risky by negative_risky_sig$PHQ_binary
## W = 461725, p-value = 0.9766
## alternative hypothesis: true location shift is not equal to 0

    cor.test(negative_risky_sig$percent_risky,negative_risky_sig$PHQ_score,method="spearman")

## Warning in cor.test.default(negative_risky_sig$percent_risky,
## negative_risky_sig$PHQ_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  negative_risky_sig$percent_risky and negative_risky_sig$PHQ_score
## S = 1142660045, p-value = 0.1319
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.03437554

  # Combine results
combined_risky <- bind_rows(positive_risky, negative_risky)

  # Plotting
  ggplot(combined_risky, aes(x = Future, y = mean_risky,  color = PHQ_binary, group = PHQ_binary)) +
    geom_point(position = position_dodge(0.7), width = 0.7) +
    geom_errorbar(aes(ymin = mean_risky - sem_risky, ymax = mean_risky + sem_risky),
                  position = position_dodge(0.7), width = 0.2) + 
    geom_hline(yintercept = 50, linetype = "dotted", color = "black") +  # Add a dotted line at y=0

    scale_color_manual(values = c("PHQ 7+" = "red", "PHQ <7" = "pink")) +
    labs(
      title = "% Risky Choices based on Future (Play 2) N=1922",
      x = "Future",
      y = "% Risky Choices"
    ) +
    ylim(40,70) +
    theme_minimal() +
    poster_theme 

## Warning in geom_point(position = position_dodge(0.7), width = 0.7): Ignoring
## unknown parameters: `width`

  # PHQ Bin
  
 positive_risky <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & !is.na(PHQ_bin)) %>%
  group_by(userKey,PHQ_bin) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
    group_by(PHQ_bin) %>%
    summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Positive")
    
    
  
  negative_risky <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & !is.na(PHQ_bin)) %>%
  group_by(userKey,PHQ_bin) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  group_by(PHQ_bin) %>%
  summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Negative")


  # Combine results
combined_risky <- bind_rows(positive_risky, negative_risky)
combined_risky <- na.omit(combined_risky)

  # Plotting
  ggplot(combined_risky, aes(x = Future, y = mean_risky,  color = PHQ_bin, group = PHQ_bin)) +
    geom_point(position = position_dodge(0.7), width = 0.7) +
    geom_errorbar(aes(ymin = mean_risky - sem_risky, ymax = mean_risky + sem_risky),
                  position = position_dodge(0.7), width = 0.2) + 
    geom_hline(yintercept = 50, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
      scale_color_manual(values = c("0-4" = "#FFC0CB", "5-9" = "#CC9AA2", "10-14" = "#99737A", "15-24" = "#664D51")) +
    labs(
      title = "% Risky Choices based on Future (Play 2) N=1922",
      x = "Future",
      y = "% Risky Choices"
    ) +
        ylim(40,70) +
    theme_minimal() +
    poster_theme 

## Warning in geom_point(position = position_dodge(0.7), width = 0.7): Ignoring
## unknown parameters: `width`

Risky With Current and Future Informaiton

current_pos_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Gain") %>%
  group_by(userKey) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  mutate(Future = "Positive",
         Current = "Gain")

current_pos_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Gain") %>%
  group_by(userKey) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  mutate(Future = "Negative",
         Current = "Gain")
  
current_neg_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2, Trial == "Loss") %>%
  group_by(userKey) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  mutate(Future = "Positive",
         Current = "Loss")

current_neg_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1., Trial == "Loss") %>%
  group_by(userKey) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  mutate(Future = "Negative",
         Current = "Loss")
  
  # Combine results
combined_risky <- bind_rows(current_pos_future_pos, current_pos_future_neg, current_neg_future_pos,current_neg_future_neg)

summary_combined <- combined_risky %>%
  group_by(Current, Future) %>%
  summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()))

## `summarise()` has grouped output by 'Current'. You can override using the
## `.groups` argument.

  # Plotting
  ggplot(summary_combined, aes(x = Future, y = mean_risky, color = Current)) +
    geom_point(position = position_dodge(0.7), width = 0.7) +
    geom_errorbar(aes(ymin = mean_risky - sem_risky, ymax = mean_risky + sem_risky),
                  position = position_dodge(0.7), width = 0.2) + 
    geom_hline(yintercept = 50, linetype = "dotted", color = "black") +  # Add a dotted line at y=0

    scale_color_manual(values = c("Gain" = "#FFD700", "Loss" = "#FF6347")) +
    labs(
      title = "% Risky Choices based on Current and Future (Play 2)",
      x = "Future",
      y = "% Risky Choices"
    ) +
    ylim(40,70) + 
    theme_minimal() +
    poster_theme 

## Warning in geom_point(position = position_dodge(0.7), width = 0.7): Ignoring
## unknown parameters: `width`

  # people choose to take risks more when the future is positive, and they are currently in a gain domain
  # people choose to take the risk the least when the future is negative and they are currently in the loss domain

Split by GAD and PHQ

#GAD
current_pos_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Gain") %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
    group_by(GAD_binary) %>%
    summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Positive",
         Current = "Gain")

#t.test
current_pos_future_pos_sig <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Gain") %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
            GAD_score = first(gad7_total),
              .groups = 'drop') 
  wilcox.test(current_pos_future_pos_sig$percent_risky~current_pos_future_pos_sig$GAD_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_pos_future_pos_sig$percent_risky by current_pos_future_pos_sig$GAD_binary
## W = 452158, p-value = 0.9014
## alternative hypothesis: true location shift is not equal to 0

  cor.test(current_pos_future_pos_sig$percent_risky,current_pos_future_pos_sig$GAD_score,method="spearman")

## Warning in cor.test.default(current_pos_future_pos_sig$percent_risky,
## current_pos_future_pos_sig$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_pos_future_pos_sig$percent_risky and current_pos_future_pos_sig$GAD_score
## S = 1148218444, p-value = 0.8797
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## 0.003468689

current_pos_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Gain") %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
    group_by(GAD_binary) %>%
    summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Negative",
         Current = "Gain")
# t.test
current_pos_future_neg_sig <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Gain") %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
            GAD_score = first(gad7_total),
              .groups = 'drop') 
  wilcox.test(current_pos_future_neg_sig$percent_risky~current_pos_future_neg_sig$GAD_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_pos_future_neg_sig$percent_risky by current_pos_future_neg_sig$GAD_binary
## W = 442598, p-value = 0.3528
## alternative hypothesis: true location shift is not equal to 0

    cor.test(current_pos_future_neg_sig$percent_risky,current_pos_future_neg_sig$GAD_score,method="spearman")

## Warning in cor.test.default(current_pos_future_neg_sig$percent_risky,
## current_pos_future_neg_sig$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_pos_future_neg_sig$percent_risky and current_pos_future_neg_sig$GAD_score
## S = 1121708464, p-value = 0.2481
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.02647653

current_neg_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2, Trial == "Loss") %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  group_by(GAD_binary) %>%
  summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Positive",
         Current = "Loss")

#t.test
current_neg_future_pos_sig <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Loss") %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
            GAD_score = first(gad7_total),
              .groups = 'drop')  
  wilcox.test(current_neg_future_pos_sig$percent_risky~current_neg_future_pos_sig$GAD_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_neg_future_pos_sig$percent_risky by current_neg_future_pos_sig$GAD_binary
## W = 459922, p-value = 0.5954
## alternative hypothesis: true location shift is not equal to 0

    cor.test(current_neg_future_pos_sig$percent_risky,current_neg_future_pos_sig$GAD_score,method="spearman")

## Warning in cor.test.default(current_neg_future_pos_sig$percent_risky,
## current_neg_future_pos_sig$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_neg_future_pos_sig$percent_risky and current_neg_future_pos_sig$GAD_score
## S = 1145235532, p-value = 0.7916
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## 0.006057539

current_neg_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1, Trial == "Loss") %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              .groups = 'drop') %>%
  group_by(GAD_binary) %>%
  summarise(mean_risky = mean(percent_risky,na.rm=TRUE),
            sem_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
            .groups = 'drop') %>%
  mutate(Future = "Negative",
         Current = "Loss")

#t.test
current_neg_future_neg_sig <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Loss") %>%
  group_by(userKey,GAD_binary) %>%
  summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
            GAD_score = first(gad7_total),
              .groups = 'drop')  
  wilcox.test(current_neg_future_neg_sig$percent_risky~current_neg_future_neg_sig$GAD_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_neg_future_neg_sig$percent_risky by current_neg_future_neg_sig$GAD_binary
## W = 455224, p-value = 0.8929
## alternative hypothesis: true location shift is not equal to 0

      cor.test(current_neg_future_neg_sig$percent_risky,current_neg_future_neg_sig$GAD_score,method="spearman")

## Warning in cor.test.default(current_neg_future_neg_sig$percent_risky,
## current_neg_future_neg_sig$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_neg_future_neg_sig$percent_risky and current_neg_future_neg_sig$GAD_score
## S = 1134139148, p-value = 0.4938
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.01568802

  # Combine results
combined_risky <- bind_rows(current_pos_future_pos, current_pos_future_neg,current_neg_future_pos,current_neg_future_neg)
combined_risky <- na.omit(combined_risky)

  # Plotting
  ggplot(combined_risky, aes(x = Future, y = mean_risky,  color = GAD_binary, group = Current)) +
    geom_point(position = position_dodge(1), width = 0.7,group="GAD_binary") +
    geom_errorbar(aes(ymin = mean_risky - sem_risky, ymax = mean_risky + sem_risky),
                  position = position_dodge(1), width = 0.2, group="GAD_binary") + 
    geom_hline(yintercept = 50, linetype = "dotted", color = "black") +  # Add a dotted line at y=0

    scale_color_manual(values = c("GAD 6+" = "#00BD9D", "GAD <6" = "#8BD7D2")) +
    labs(
      title = "% Risky Choices based on Future (Play 1) N=5393",
      x = "Future",
      y = "% Risky Choices"
    ) +
    facet_wrap(~Current) +
    ylim(40,70) + 
    theme_minimal() +
    poster_theme 

## Warning in geom_point(position = position_dodge(1), width = 0.7, group =
## "GAD_binary"): Ignoring unknown parameters: `width`

Within Samples Analysis

# Comparing risky choices when future is positive versus negative within GAD GROUPS in GAIN DOMAIN
gain_data <- happyData1 %>%
  filter(Trial == "Gain", !is.na(GAD_binary)) %>%
  group_by(userKey, GAD_binary, Future = ifelse(NextIsland == 2 & LastIsland == 2, "Positive", "Negative")) %>%
  summarise(
    total_choices = n(),
    risky_choices = sum(Choice == 2, na.rm = TRUE),
    percent_risky = (risky_choices / total_choices) * 100,
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Future, values_from = percent_risky)

gain_ANX <- gain_data %>%
  filter(GAD_binary == "GAD 6+")

gain_ANX_long <- gain_ANX %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "percent_risky"     # New column for percent risky choices
  )

wilcox.test(gain_ANX_long$percent_risky~gain_ANX_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  gain_ANX_long$percent_risky by gain_ANX_long$Future
## W = 428411, p-value = 0.02643
## alternative hypothesis: true location shift is not equal to 0

gain_lowANX <- gain_data %>%
  filter(GAD_binary == "GAD <6")

gain_lowANX_long <- gain_lowANX %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "percent_risky"     # New column for percent risky choices
  )


wilcox.test(gain_lowANX_long$percent_risky~gain_lowANX_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  gain_lowANX_long$percent_risky by gain_lowANX_long$Future
## W = 418014, p-value = 0.004214
## alternative hypothesis: true location shift is not equal to 0

# IN LOSS DOMAIN
loss_data <- happyData1 %>%
  filter(Trial == "Loss", !is.na(GAD_binary)) %>%
  group_by(userKey, GAD_binary, Future = ifelse(NextIsland == 2 & LastIsland == 2, "Positive", "Negative")) %>%
  summarise(
    total_choices = n(),
    risky_choices = sum(Choice == 2, na.rm = TRUE),
    percent_risky = (risky_choices / total_choices) * 100,
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Future, values_from = percent_risky)

loss_ANX <- loss_data %>%
  filter(GAD_binary == "GAD 6+")

loss_ANX_long <- loss_ANX %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "percent_risky"     # New column for percent risky choices
  )

wilcox.test(loss_ANX_long$percent_risky~loss_ANX_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  loss_ANX_long$percent_risky by loss_ANX_long$Future
## W = 437475, p-value = 0.1428
## alternative hypothesis: true location shift is not equal to 0

loss_lowANX <- loss_data %>%
  filter(GAD_binary == "GAD <6")

loss_lowANX_long <- loss_lowANX %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "percent_risky"     # New column for percent risky choices
  )


wilcox.test(loss_lowANX_long$percent_risky~loss_lowANX_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  loss_lowANX_long$percent_risky by loss_lowANX_long$Future
## W = 425910, p-value = 0.02769
## alternative hypothesis: true location shift is not equal to 0

# comparing % risky choices when future is positive in Gain vs Loss Domains

positive_future_data <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2, !is.na(GAD_binary)) %>%  # Filter for Positive future
  group_by(userKey, GAD_binary, Domain = ifelse(Trial == "Gain", "Gain", "Loss")) %>%
  summarise(
    total_choices = n(),
    risky_choices = sum(Choice == 2, na.rm = TRUE),
    percent_risky = (risky_choices / total_choices) * 100,
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Domain, values_from = percent_risky)

positive_future_long <- positive_future_data %>%
  pivot_longer(
    cols = c(Gain, Loss),  # Columns to pivot
    names_to = "Domain",   # New column for domain condition (Gain vs Loss)
    values_to = "percent_risky"  # New column for percent risky choices
  ) %>%
  drop_na(percent_risky)  # Remove any rows with NA values in percent_risky

positive_future_anx <- positive_future_long %>%
  filter(GAD_binary == "GAD 6+")

wilcox.test(percent_risky ~ Domain, data = positive_future_anx)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  percent_risky by Domain
## W = 436490, p-value = 0.1184
## alternative hypothesis: true location shift is not equal to 0

# Repeat for GAD <5 group
positive_future_non_anx <- positive_future_long %>%
  filter(GAD_binary == "GAD <6")

wilcox.test(percent_risky ~ Domain,  data = positive_future_non_anx)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  percent_risky by Domain
## W = 423830, p-value = 0.01655
## alternative hypothesis: true location shift is not equal to 0

## for negative future
negative_future_data <- happyData1 %>%
  filter(NextIsland == 1, LastIsland == 1, !is.na(GAD_binary)) %>%  # Filter for Positive future
  group_by(userKey, GAD_binary, Domain = ifelse(Trial == "Gain", "Gain", "Loss")) %>%
  summarise(
    total_choices = n(),
    risky_choices = sum(Choice == 2, na.rm = TRUE),
    percent_risky = (risky_choices / total_choices) * 100,
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Domain, values_from = percent_risky)

negative_future_long <- negative_future_data %>%
  pivot_longer(
    cols = c(Gain, Loss),  # Columns to pivot
    names_to = "Domain",   # New column for domain condition (Gain vs Loss)
    values_to = "percent_risky"  # New column for percent risky choices
  ) %>%
  drop_na(percent_risky)  # Remove any rows with NA values in percent_risky

negative_future_anx <- negative_future_long %>%
  filter(GAD_binary == "GAD 6+")

wilcox.test(percent_risky ~ Domain, data = negative_future_anx)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  percent_risky by Domain
## W = 437284, p-value = 0.1352
## alternative hypothesis: true location shift is not equal to 0

# Repeat for GAD <5 group
negative_future_non_anx <- negative_future_long %>%
  filter(GAD_binary == "GAD <6")

wilcox.test(percent_risky ~ Domain,  data = negative_future_non_anx)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  percent_risky by Domain
## W = 421286, p-value = 0.009024
## alternative hypothesis: true location shift is not equal to 0

Future and Current Information Happiness

current_pos_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Gain") %>%
  group_by(userKey) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE)) %>%
  mutate(Future = "Positive",
         Current = "Gain")

current_pos_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Gain") %>%
  group_by(userKey) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE)) %>%
  mutate(Future = "Negative",
         Current = "Gain")
  
current_neg_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2, Trial == "Loss") %>%
  group_by(userKey) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE)) %>%
  mutate(Future = "Positive",
         Current = "Loss")

current_neg_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1., Trial == "Loss") %>%
  group_by(userKey) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE)) %>%
  mutate(Future = "Negative",
         Current = "Loss")
  
  # Combine results
combined_happy <- bind_rows(current_pos_future_pos, current_pos_future_neg, current_neg_future_pos,current_neg_future_neg)


 ggplot(combined_happy, aes(x = Future, y = overall_mean_happiness, fill = Current)) +
    geom_bar(stat = "identity", position = position_dodge(), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = overall_mean_happiness - sem_happiness, ymax = overall_mean_happiness + sem_happiness),
                  position = position_dodge(0.7), width = 0.2) +
    geom_point(aes(y = overall_mean_pred_happiness), color = "lightblue", shape = 8, size = 3, position=position_dodge(0.7)) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0

    scale_fill_manual(values = c("Gain" = "#FFD700", "Loss" = "#FF6347")) +
    labs(
      title = "z-scored Happiness based on Future and Current",
      x = "Future",
      y = "z-scored Happiness"
    ) +
    theme_minimal() +
    poster_theme 

  # people choose to take risks more when the future is positive, and they are currently in a gain domain
  # people choose to take the risk the least when the future is negative and they are currently in the loss domain

combined_happy <- combined_happy %>%
  mutate(residuals = overall_mean_pred_happiness - overall_mean_happiness)

# Plot residuals with error bars
ggplot(combined_happy, aes(x = Future, y = residuals, color = Current)) +
  geom_point(aes(y = residuals), size = 3, position = position_dodge(0.7)) +  # Residual dots
  geom_errorbar(aes(ymin = residuals - sem_happiness, ymax = residuals + sem_happiness),
                position = position_dodge(0.7), width = 0.2) +  # Error bars for residuals
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
  
  scale_color_manual(values = c("Gain" = "#FFD700", "Loss" = "#FF6347")) +
  labs(
    title = "Residuals of Predicted vs Actual z-scored Happiness",
    x = "Future",
    y = "Residuals (Predicted - Actual Happiness)"
  ) +
  theme_minimal() +
  ylim(-0.05, 0.05) +
  poster_theme

Happiness Based on Future and Current - Split by GAD and PHQ

current_pos_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Gain") %>%
  group_by(userKey, GAD_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Positive",
         Current = "Gain")

#t.test
current_pos_future_pos_sig <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Gain") %>%
  group_by(userKey, GAD_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            GAD_score=first(gad7_total),
            .groups = 'drop'
            )
  wilcox.test(current_pos_future_pos_sig$mean_happiness~current_pos_future_pos_sig$GAD_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_pos_future_pos_sig$mean_happiness by current_pos_future_pos_sig$GAD_binary
## W = 241037, p-value = 0.5615
## alternative hypothesis: true location shift is not equal to 0

      cor.test(current_pos_future_pos_sig$mean_happiness,current_pos_future_pos_sig$GAD_score,method="spearman")

## Warning in cor.test.default(current_pos_future_pos_sig$mean_happiness,
## current_pos_future_pos_sig$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_pos_future_pos_sig$mean_happiness and current_pos_future_pos_sig$GAD_score
## S = 460935716, p-value = 0.8937
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.00357028

current_pos_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Gain") %>%
  group_by(userKey, GAD_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Negative",
         Current = "Gain")

#t.test
current_pos_future_neg_sig <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Gain") %>%
  group_by(userKey, GAD_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            GAD_score=first(gad7_total),
            .groups = 'drop')
  wilcox.test(current_pos_future_neg_sig$mean_happiness~current_pos_future_neg_sig$GAD_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_pos_future_neg_sig$mean_happiness by current_pos_future_neg_sig$GAD_binary
## W = 216531, p-value = 0.08024
## alternative hypothesis: true location shift is not equal to 0

        cor.test(current_pos_future_neg_sig$mean_happiness,current_pos_future_neg_sig$GAD_score,method="spearman")

## Warning in cor.test.default(current_pos_future_neg_sig$mean_happiness,
## current_pos_future_neg_sig$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_pos_future_neg_sig$mean_happiness and current_pos_future_neg_sig$GAD_score
## S = 388849369, p-value = 0.02698
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.06011029

current_neg_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2, Trial == "Loss") %>%
  group_by(userKey, GAD_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Positive",
         Current = "Loss")

#t.test
current_neg_future_pos_sig <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Loss") %>%
  group_by(userKey, GAD_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            GAD_score=first(gad7_total),
            .groups = 'drop')
  wilcox.test(current_neg_future_pos_sig$mean_happiness~current_neg_future_pos_sig$GAD_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_neg_future_pos_sig$mean_happiness by current_neg_future_pos_sig$GAD_binary
## W = 237840, p-value = 0.1753
## alternative hypothesis: true location shift is not equal to 0

        cor.test(current_neg_future_pos_sig$mean_happiness,current_neg_future_pos_sig$GAD_score,method="spearman")

## Warning in cor.test.default(current_neg_future_pos_sig$mean_happiness,
## current_neg_future_pos_sig$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_neg_future_pos_sig$mean_happiness and current_neg_future_pos_sig$GAD_score
## S = 433772879, p-value = 0.04148
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.05547472

current_neg_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1., Trial == "Loss") %>%
  group_by(userKey, GAD_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Negative",
         Current = "Loss")

#t.test
current_neg_future_neg_sig <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Loss") %>%
  group_by(userKey, GAD_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            GAD_score=first(gad7_total),
            .groups = 'drop')
  wilcox.test(current_neg_future_neg_sig$mean_happiness~current_neg_future_neg_sig$GAD_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_neg_future_neg_sig$mean_happiness by current_neg_future_neg_sig$GAD_binary
## W = 269534, p-value = 0.5629
## alternative hypothesis: true location shift is not equal to 0

          cor.test(current_neg_future_neg_sig$mean_happiness,current_neg_future_neg_sig$GAD_score,method="spearman")

## Warning in cor.test.default(current_neg_future_neg_sig$mean_happiness,
## current_neg_future_neg_sig$GAD_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_neg_future_neg_sig$mean_happiness and current_neg_future_neg_sig$GAD_score
## S = 516169948, p-value = 0.8978
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##          rho 
## -0.003367538

  # Combine results
combined_happy <- bind_rows(current_pos_future_pos, current_pos_future_neg, current_neg_future_pos,current_neg_future_neg)
combined_happy <- na.omit(combined_happy)

 ggplot(combined_happy, aes(x = Future, y = overall_mean_happiness, fill = GAD_binary, group=GAD_binary)) +
    geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = overall_mean_happiness - sem_happiness, ymax = overall_mean_happiness + sem_happiness),
                  position = position_dodge(0.7), width = 0.2) +
    geom_point(aes(y = overall_mean_pred_happiness), color = "lightblue", shape = 8, size = 3, position=position_dodge(0.7)) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    scale_fill_manual(values = c("GAD 6+" = "#00BD9D", "GAD <6" = "#8BD7D2")) +
    labs(
      title = "z-scored Happiness based on Future and Current",
      x = "Future",
      y = "z-scored Happiness"
    ) +
   facet_wrap(~Current) +
    theme_minimal() +
    poster_theme 

 ## GAD BINS - fill= Current, wrap~GAD_bin
 
 current_pos_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Gain") %>%
  group_by(userKey, GAD_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Positive",
         Current = "Gain")
 
 current_pos_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Gain") %>%
  group_by(userKey, GAD_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Negative",
         Current = "Gain")
 
 current_neg_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2, Trial == "Loss") %>%
  group_by(userKey, GAD_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Positive",
         Current = "Loss")
 
 current_neg_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1., Trial == "Loss") %>%
  group_by(userKey, GAD_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(GAD_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Negative",
         Current = "Loss")
 
 # Combine results
combined_happy <- bind_rows(current_pos_future_pos, current_pos_future_neg, current_neg_future_pos,current_neg_future_neg)
combined_happy <- na.omit(combined_happy)

 ggplot(combined_happy, aes(x = Future, y = overall_mean_happiness, fill = GAD_bin, group=GAD_bin)) +
    geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = overall_mean_happiness - sem_happiness, ymax = overall_mean_happiness + sem_happiness),
                  position = position_dodge(0.7), width = 0.2) +
    geom_point(aes(y = overall_mean_pred_happiness), color = "lightblue", shape = 8, size = 3, position=position_dodge(0.7)) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    scale_fill_manual(values = c("0-4" = "#8BD7D2", "5-9" = "#6FACA8", "10-14" = "#53817E", "15-21" = "#385654")) +
    labs(
      title = "z-scored Happiness based on Future and Current",
      x = "Future",
      y = "z-scored Happiness"
    ) +
   facet_grid(~Current) +
    theme_minimal() +
    poster_theme 

 ### PHQ
 current_pos_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Gain") %>%
  group_by(userKey, PHQ_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Positive",
         Current = "Gain")

#t.test
current_pos_future_pos_sig <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Gain") %>%
  group_by(userKey, PHQ_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            PHQ_score=first(phq8_total),
            .groups = 'drop'
            )
  wilcox.test(current_pos_future_pos_sig$mean_happiness~current_pos_future_pos_sig$PHQ_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_pos_future_pos_sig$mean_happiness by current_pos_future_pos_sig$PHQ_binary
## W = 251020, p-value = 0.9248
## alternative hypothesis: true location shift is not equal to 0

      cor.test(current_pos_future_pos_sig$mean_happiness,current_pos_future_pos_sig$PHQ_score,method="spearman")

## Warning in cor.test.default(current_pos_future_pos_sig$mean_happiness,
## current_pos_future_pos_sig$PHQ_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_pos_future_pos_sig$mean_happiness and current_pos_future_pos_sig$PHQ_score
## S = 482700289, p-value = 0.45
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.02009182

current_pos_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Gain") %>%
  group_by(userKey, PHQ_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Negative",
         Current = "Gain")

#t.test
current_pos_future_neg_sig <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Gain") %>%
  group_by(userKey, PHQ_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            PHQ_score=first(phq8_total),
            .groups = 'drop')
  wilcox.test(current_pos_future_neg_sig$mean_happiness~current_pos_future_neg_sig$PHQ_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_pos_future_neg_sig$mean_happiness by current_pos_future_neg_sig$PHQ_binary
## W = 217988, p-value = 0.05665
## alternative hypothesis: true location shift is not equal to 0

        cor.test(current_pos_future_neg_sig$mean_happiness,current_pos_future_neg_sig$PHQ_score,method="spearman")

## Warning in cor.test.default(current_pos_future_neg_sig$mean_happiness,
## current_pos_future_neg_sig$PHQ_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_pos_future_neg_sig$mean_happiness and current_pos_future_neg_sig$PHQ_score
## S = 402170735, p-value = 0.09735
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.04494006

current_neg_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2, Trial == "Loss") %>%
  group_by(userKey, PHQ_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Positive",
         Current = "Loss")

#t.test
current_neg_future_pos_sig <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Loss") %>%
  group_by(userKey, PHQ_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            PHQ_score=first(phq8_total),
            .groups = 'drop')
  wilcox.test(current_neg_future_pos_sig$mean_happiness~current_neg_future_pos_sig$PHQ_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_neg_future_pos_sig$mean_happiness by current_neg_future_pos_sig$PHQ_binary
## W = 244015, p-value = 0.1229
## alternative hypothesis: true location shift is not equal to 0

        cor.test(current_neg_future_pos_sig$mean_happiness,current_neg_future_pos_sig$PHQ_score,method="spearman")

## Warning in cor.test.default(current_neg_future_pos_sig$mean_happiness,
## current_neg_future_pos_sig$PHQ_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_neg_future_pos_sig$mean_happiness and current_neg_future_pos_sig$PHQ_score
## S = 447264725, p-value = 0.04159
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.05515965

current_neg_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1., Trial == "Loss") %>%
  group_by(userKey, PHQ_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_binary) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Negative",
         Current = "Loss")

#t.test
current_neg_future_neg_sig <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Loss") %>%
  group_by(userKey, PHQ_binary) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            PHQ_score=first(phq8_total),
            .groups = 'drop')
  wilcox.test(current_neg_future_neg_sig$mean_happiness~current_neg_future_neg_sig$PHQ_binary)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  current_neg_future_neg_sig$mean_happiness by current_neg_future_neg_sig$PHQ_binary
## W = 271378, p-value = 0.9305
## alternative hypothesis: true location shift is not equal to 0

          cor.test(current_neg_future_neg_sig$mean_happiness,current_neg_future_neg_sig$PHQ_score,method="spearman")

## Warning in cor.test.default(current_neg_future_neg_sig$mean_happiness,
## current_neg_future_neg_sig$PHQ_score, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  current_neg_future_neg_sig$mean_happiness and current_neg_future_neg_sig$PHQ_score
## S = 530400031, p-value = 0.932
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## 0.002227427

  # Combine results
combined_happy <- bind_rows(current_pos_future_pos, current_pos_future_neg, current_neg_future_pos,current_neg_future_neg)
combined_happy <- na.omit(combined_happy)

 ggplot(combined_happy, aes(x = Future, y = overall_mean_happiness, fill = PHQ_binary, group=PHQ_binary)) +
    geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = overall_mean_happiness - sem_happiness, ymax = overall_mean_happiness + sem_happiness),
                  position = position_dodge(0.7), width = 0.2) +
    geom_point(aes(y = overall_mean_pred_happiness), color = "lightblue", shape = 8, size = 3, position=position_dodge(0.7)) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    scale_fill_manual(values = c("PHQ 7+" = "red", "PHQ <7" = "pink")) +
    labs(
      title = "z-scored Happiness based on Future and Current",
      x = "Future",
      y = "z-scored Happiness"
    ) +
   facet_wrap(~Current) +
    theme_minimal() +
    poster_theme 

 ## GAD BINS - fill= Current, wrap~PHQ_bin
 
 current_pos_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2 & Trial == "Gain") %>%
  group_by(userKey, PHQ_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Positive",
         Current = "Gain")
 
 current_pos_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1 & Trial == "Gain") %>%
  group_by(userKey, PHQ_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Negative",
         Current = "Gain")
 
 current_neg_future_pos <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2, Trial == "Loss") %>%
  group_by(userKey, PHQ_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Positive",
         Current = "Loss")
 
 current_neg_future_neg <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1., Trial == "Loss") %>%
  group_by(userKey, PHQ_bin) %>%
  summarise(mean_happiness = mean(zHappy, na.rm = TRUE),
            mean_pred_happiness = mean(zHappyPred, na.rm = TRUE), 
            .groups = 'drop') %>%
  group_by(PHQ_bin) %>%
  summarise(overall_mean_happiness = mean(mean_happiness, na.rm = TRUE),
            sem_happiness = sd(mean_happiness, na.rm = TRUE) / sqrt(n()),
            overall_mean_pred_happiness = mean(mean_pred_happiness, na.rm = TRUE),
            .groups = 'drop')  %>%
  mutate(Future = "Negative",
         Current = "Loss")
 
 # Combine results
combined_happy <- bind_rows(current_pos_future_pos, current_pos_future_neg, current_neg_future_pos,current_neg_future_neg)
combined_happy <- na.omit(combined_happy)

 ggplot(combined_happy, aes(x = Future, y = overall_mean_happiness, fill = PHQ_bin, group=PHQ_bin)) +
    geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = overall_mean_happiness - sem_happiness, ymax = overall_mean_happiness + sem_happiness),
                  position = position_dodge(0.7), width = 0.2) +
    geom_point(aes(y = overall_mean_pred_happiness), color = "lightblue", shape = 8, size = 3, position=position_dodge(0.7)) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
      scale_fill_manual(values = c("0-4" = "#FFC0CB", "5-9" = "#CC9AA2", "10-14" = "#99737A", "15-24" = "#664D51")) +
    labs(
      title = "z-scored Happiness based on Future and Current",
      x = "Future",
      y = "z-scored Happiness"
    ) +
   facet_grid(~Current) +
    theme_minimal() +
    poster_theme 

Within Groups Analysis

# Comparing happiness when future is positive versus negative within GAD GROUPS in GAIN DOMAIN
gain_data <- happyData1 %>%
  filter(Trial == "Gain", !is.na(GAD_binary)) %>%
  group_by(userKey, GAD_binary, Future = ifelse(NextIsland == 2 & LastIsland == 2, "Positive", "Negative")) %>%
  summarise(
    mean_happiness = mean(zHappy, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Future, values_from = mean_happiness)

gain_ANX <- gain_data %>%
  filter(GAD_binary == "GAD 6+")

gain_ANX_long <- gain_ANX %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "mean_happiness"     # New column for percent risky choices
  )


wilcox.test(gain_ANX_long$mean_happiness~gain_ANX_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  gain_ANX_long$mean_happiness by gain_ANX_long$Future
## W = 363642, p-value = 0.06275
## alternative hypothesis: true location shift is not equal to 0

gain_lowANX <- gain_data %>%
  filter(GAD_binary == "GAD <6")

gain_lowANX_long <- gain_lowANX %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "mean_happiness"     # New column for percent risky choices
  )

wilcox.test(gain_lowANX_long$mean_happiness~gain_lowANX_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  gain_lowANX_long$mean_happiness by gain_lowANX_long$Future
## W = 334422, p-value = 0.1985
## alternative hypothesis: true location shift is not equal to 0

#loss
loss_data <- happyData1 %>%
  filter(Trial == "Loss", !is.na(GAD_binary)) %>%
  group_by(userKey, GAD_binary, Future = ifelse(NextIsland == 2 & LastIsland == 2, "Positive", "Negative")) %>%
  summarise(
    mean_happiness = mean(zHappy, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Future, values_from = mean_happiness)

loss_ANX <- loss_data %>%
  filter(GAD_binary == "GAD 6+")

loss_ANX_long <- loss_ANX %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "mean_happiness"     # New column for percent risky choices
  )


wilcox.test(loss_ANX_long$mean_happiness~loss_ANX_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  loss_ANX_long$mean_happiness by loss_ANX_long$Future
## W = 343600, p-value = 0.007168
## alternative hypothesis: true location shift is not equal to 0

loss_lowANX <- loss_data %>%
  filter(GAD_binary == "GAD <6")

loss_lowANX_long <- loss_lowANX %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "mean_happiness"     # New column for percent risky choices
  )

wilcox.test(loss_lowANX_long$mean_happiness~loss_lowANX_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  loss_lowANX_long$mean_happiness by loss_lowANX_long$Future
## W = 337338, p-value = 0.1814
## alternative hypothesis: true location shift is not equal to 0

# Comparing happiness when future is Positive in Gain vs Loss Domain

positive_future_data <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2, !is.na(GAD_binary)) %>%  # Filter for Positive future
  group_by(userKey, GAD_binary, Domain = ifelse(Trial == "Gain", "Gain", "Loss")) %>%
  summarise(
   mean_happiness = mean(zHappy, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Domain, values_from = mean_happiness)

positive_future_long <- positive_future_data %>%
  pivot_longer(
    cols = c(Gain, Loss),  # Columns to pivot
    names_to = "Domain",   # New column for domain condition (Gain vs Loss)
    values_to = "mean_happiness"  # New column for percent risky choices
  ) %>%
  drop_na(mean_happiness)  # Remove any rows with NA values in percent_risky

positive_future_anx <- positive_future_long %>%
  filter(GAD_binary == "GAD 6+")

wilcox.test(mean_happiness ~ Domain, data = positive_future_anx)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mean_happiness by Domain
## W = 300222, p-value = 6.44e-15
## alternative hypothesis: true location shift is not equal to 0

# Repeat for GAD <5 group
positive_future_non_anx <- positive_future_long %>%
  filter(GAD_binary == "GAD <6")

wilcox.test(mean_happiness ~ Domain,  data = positive_future_non_anx)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mean_happiness by Domain
## W = 272849, p-value = 1.209e-08
## alternative hypothesis: true location shift is not equal to 0

# When future is negaative in gain vs loss domain
negative_future_data <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1, !is.na(GAD_binary)) %>%  # Filter for Positive future
  group_by(userKey, GAD_binary, Domain = ifelse(Trial == "Gain", "Gain", "Loss")) %>%
  summarise(
   mean_happiness = mean(zHappy, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Domain, values_from = mean_happiness)

negative_future_long <- negative_future_data %>%
  pivot_longer(
    cols = c(Gain, Loss),  # Columns to pivot
    names_to = "Domain",   # New column for domain condition (Gain vs Loss)
    values_to = "mean_happiness"  # New column for percent risky choices
  ) %>%
  drop_na(mean_happiness)  # Remove any rows with NA values in percent_risky

negative_future_anx <- negative_future_long %>%
  filter(GAD_binary == "GAD 6+")

wilcox.test(mean_happiness ~ Domain, data = negative_future_anx)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mean_happiness by Domain
## W = 297176, p-value = 7.034e-11
## alternative hypothesis: true location shift is not equal to 0

# Repeat for GAD <5 group
negative_future_non_anx <- negative_future_long %>%
  filter(GAD_binary == "GAD <6")

wilcox.test(mean_happiness ~ Domain,  data = negative_future_non_anx)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mean_happiness by Domain
## W = 279234, p-value = 7.19e-06
## alternative hypothesis: true location shift is not equal to 0

### PHQ 

# Comparing happiness when future is positive versus negative within GAD GROUPS in GAIN DOMAIN
gain_data <- happyData1 %>%
  filter(Trial == "Gain", !is.na(PHQ_binary)) %>%
  group_by(userKey, PHQ_binary, Future = ifelse(NextIsland == 2 & LastIsland == 2, "Positive", "Negative")) %>%
  summarise(
    mean_happiness = mean(zHappy, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Future, values_from = mean_happiness)

gain_DEP <- gain_data %>%
  filter(PHQ_binary == "PHQ 7+")

gain_DEP_long <- gain_DEP %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "mean_happiness"     # New column for percent risky choices
  )


wilcox.test(gain_DEP_long$mean_happiness~gain_DEP_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  gain_DEP_long$mean_happiness by gain_DEP_long$Future
## W = 385571, p-value = 0.02684
## alternative hypothesis: true location shift is not equal to 0

gain_lowDEP <- gain_data %>%
  filter(PHQ_binary == "PHQ <7")

gain_lowDEP_long <- gain_lowDEP %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "mean_happiness"     # New column for percent risky choices
  )

wilcox.test(gain_lowDEP_long$mean_happiness~gain_lowDEP_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  gain_lowDEP_long$mean_happiness by gain_lowDEP_long$Future
## W = 327256, p-value = 0.3254
## alternative hypothesis: true location shift is not equal to 0

#loss
loss_data <- happyData1 %>%
  filter(Trial == "Loss", !is.na(PHQ_binary)) %>%
  group_by(userKey, PHQ_binary, Future = ifelse(NextIsland == 2 & LastIsland == 2, "Positive", "Negative")) %>%
  summarise(
    mean_happiness = mean(zHappy, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Future, values_from = mean_happiness)

loss_DEP <- loss_data %>%
  filter(PHQ_binary == "PHQ 7+")

loss_DEP_long <- loss_DEP %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "mean_happiness"     # New column for percent risky choices
  )


wilcox.test(loss_DEP_long$mean_happiness~loss_DEP_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  loss_DEP_long$mean_happiness by loss_DEP_long$Future
## W = 374468, p-value = 0.002939
## alternative hypothesis: true location shift is not equal to 0

loss_lowDEP <- loss_data %>%
  filter(PHQ_binary == "PHQ <7")

loss_lowDEP_long <- loss_lowDEP %>%
  pivot_longer(
    cols = c(Positive, Negative),  # Columns to pivot
    names_to = "Future",            # New column for future condition
    values_to = "mean_happiness"     # New column for percent risky choices
  )

wilcox.test(loss_lowDEP_long$mean_happiness~loss_lowDEP_long$Future)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  loss_lowDEP_long$mean_happiness by loss_lowDEP_long$Future
## W = 320092, p-value = 0.3267
## alternative hypothesis: true location shift is not equal to 0

# Comparing happiness when future is Positive in Gain vs Loss Domain

positive_future_data <- happyData1 %>%
  filter(NextIsland == 2 & LastIsland == 2, !is.na(PHQ_binary)) %>%  # Filter for Positive future
  group_by(userKey, PHQ_binary, Domain = ifelse(Trial == "Gain", "Gain", "Loss")) %>%
  summarise(
   mean_happiness = mean(zHappy, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Domain, values_from = mean_happiness)

positive_future_long <- positive_future_data %>%
  pivot_longer(
    cols = c(Gain, Loss),  # Columns to pivot
    names_to = "Domain",   # New column for domain condition (Gain vs Loss)
    values_to = "mean_happiness"  # New column for percent risky choices
  ) %>%
  drop_na(mean_happiness)  # Remove any rows with NA values in percent_risky

positive_future_dep <- positive_future_long %>%
  filter(PHQ_binary == "PHQ 7+")

wilcox.test(mean_happiness ~ Domain, data = positive_future_dep)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mean_happiness by Domain
## W = 315198, p-value = 4.449e-14
## alternative hypothesis: true location shift is not equal to 0

# Repeat for GAD <5 group
positive_future_non_dep <- positive_future_long %>%
  filter(PHQ_binary == "PHQ <7")

wilcox.test(mean_happiness ~ Domain,  data = positive_future_non_dep)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mean_happiness by Domain
## W = 270604, p-value = 1.612e-09
## alternative hypothesis: true location shift is not equal to 0

# When future is negaative in gain vs loss domain
negative_future_data <- happyData1 %>%
  filter(NextIsland == 1 & LastIsland == 1, !is.na(PHQ_binary)) %>%  # Filter for Positive future
  group_by(userKey, PHQ_binary, Domain = ifelse(Trial == "Gain", "Gain", "Loss")) %>%
  summarise(
   mean_happiness = mean(zHappy, na.rm = TRUE),
    .groups = 'drop'
  ) %>%
  pivot_wider(names_from = Domain, values_from = mean_happiness)

negative_future_long <- negative_future_data %>%
  pivot_longer(
    cols = c(Gain, Loss),  # Columns to pivot
    names_to = "Domain",   # New column for domain condition (Gain vs Loss)
    values_to = "mean_happiness"  # New column for percent risky choices
  ) %>%
  drop_na(mean_happiness)  # Remove any rows with NA values in percent_risky

negative_future_dep <- negative_future_long %>%
  filter(PHQ_binary == "PHQ 7+")

wilcox.test(mean_happiness ~ Domain, data = negative_future_dep)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mean_happiness by Domain
## W = 312608, p-value = 9.749e-11
## alternative hypothesis: true location shift is not equal to 0

# Repeat for GAD <5 group
negative_future_non_dep <- negative_future_long %>%
  filter(PHQ_binary == "PHQ <7")

wilcox.test(mean_happiness ~ Domain,  data = negative_future_non_dep)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  mean_happiness by Domain
## W = 274216, p-value = 5.017e-06
## alternative hypothesis: true location shift is not equal to 0

Happiness Based on Gamble Decisions

#z-scored happiness
choice_risky <- happyData1 %>%
  filter(Choice == 2) %>%
  group_by(userKey) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky")
  
# safe choice
choice_safe <- happyData1 %>%
  filter(Choice == 1) %>%
  group_by(userKey) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe")
  
choice_combined <- bind_rows(choice_risky, choice_safe)
summary_choice <- choice_combined %>%
    group_by(Choice) %>% 
  summarise(mean_happy = mean(mean_zhappy,na.rm=TRUE),
            mean_predhappy = mean(mean_zpredhappy,na.rm=TRUE),
            sem_happy = sd(mean_zhappy, na.rm = TRUE) / sqrt(n()),
            sem_predhappy = sd(mean_zpredhappy, na.rm = TRUE) / sqrt(n()))


ggplot(summary_choice, aes(x = Choice, y = mean_happy,  fill = Choice, group = Choice)) +
    geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = mean_happy - sem_happy, ymax = mean_happy + sem_happy),
                  position = position_dodge(1), width = 0.2) + 
      geom_point(aes(y = mean_predhappy), color = "lightblue", shape = 8, size = 3, position=position_dodge(0.7)) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    labs(
      title = "Happiness  (z-scored) Based on Choice",
      x = "Choice",
      y = "zHappy"
    ) +
    theme_minimal() +
  ylim(-0.03, 0.03) +
    poster_theme 

# residuals
# happiness based on risky vs safe choice (backfilled)
happyData1 <- happyData1 %>%
  mutate(residuals=zHappy_filled - zHappyPred_filled)

choice_risky <- happyData1 %>%
  filter(Choice == 2) %>%
  group_by(userKey) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky")
  
# safe choice
choice_safe <- happyData1 %>%
  filter(Choice == 1) %>%
  group_by(userKey) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe")
  
choice_combined <- bind_rows(choice_risky, choice_safe)
summary_choice <- choice_combined %>%
    group_by(Choice) %>% 
  summarise(mean_resids = mean(mean_residuals,na.rm=TRUE),
            sem_resids = sd(mean_residuals, na.rm = TRUE) / sqrt(n()))


ggplot(summary_choice, aes(x = Choice, y = mean_resids,  fill = Choice, grgoup=Choice)) +
    geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = mean_resids - sem_resids, ymax = mean_resids + sem_resids),
                  position = position_dodge(1), width = 0.2) + 
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    labs(
      title = "Happiness Residuals (z-scored) Based on Choice",
      x = "Choice",
      y = "zHappy Resid"
    ) +
    theme_minimal() +
  ylim(-0.03, 0.03) + 
    poster_theme 

Happiness Based on Choice and Current Domain

happyData1 <- happyData1 %>%
  mutate(residuals=zHappy_filled - zHappyPred_filled)

#z-scored happiness
choice_risky_gain <- happyData1 %>%
  filter(Choice == 2 & Trial =="Gain") %>%
  group_by(userKey) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Current = "Gain")


choice_risky_loss <- happyData1 %>%
  filter(Choice == 2 & Trial =="Loss") %>%
  group_by(userKey) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Current = "Loss")
  
# safe choice
choice_safe_gain <- happyData1 %>%
  filter(Choice == 1 & Trial=="Gain") %>%
  group_by(userKey) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
        Current = "Gain")

choice_safe_loss <- happyData1 %>%
  filter(Choice == 1 & Trial=="Loss") %>%
  group_by(userKey) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
        Current = "Loss")

  
choice_combined <- bind_rows(choice_risky_gain, choice_risky_loss,choice_safe_gain,choice_safe_loss)
summary_choice <- choice_combined %>%
    group_by(Choice,Current) %>% 
  summarise(mean_happy = mean(mean_zhappy,na.rm=TRUE),
            mean_predhappy = mean(mean_zpredhappy,na.rm=TRUE),
            sem_happy = sd(mean_zhappy, na.rm = TRUE) / sqrt(n()),
            sem_predhappy = sd(mean_zpredhappy, na.rm = TRUE) / sqrt(n()))

## `summarise()` has grouped output by 'Choice'. You can override using the
## `.groups` argument.

ggplot(summary_choice, aes(x = Choice, y = mean_happy,  fill = Current, group=Current)) +
    geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = mean_happy - sem_happy, ymax = mean_happy + sem_happy),position = position_dodge(0.7), width = 0.2) + 
      geom_point(aes(y = mean_predhappy), color = "lightblue", shape = 8, size = 3,position = position_dodge(0.7)) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    scale_fill_manual(values = c("Gain" = "#FFD700", "Loss" = "#FF6347")) +
    labs(
      title = "Happiness  (z-scored) Based on Choice",
      x = "Choice",
      y = "zHappy"
    ) +
    theme_minimal() +
    poster_theme 

# residuals

choice_risky_gain <- happyData1 %>%
  filter(Choice == 2 & Trial == "Gain") %>%
  group_by(userKey) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Current = "Gain")

choice_risky_loss <- happyData1 %>%
  filter(Choice == 2 & Trial == "Loss") %>%
  group_by(userKey) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Current = "Loss")


  
# safe choice
choice_safe_gain <- happyData1 %>%
  filter(Choice == 1 & Trial == "Gain") %>%
  group_by(userKey) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
                  Current = "Gain")

choice_safe_loss <- happyData1 %>%
  filter(Choice == 1 & Trial == "Loss") %>%
  group_by(userKey) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
                  Current = "Loss")

  

choice_combined <- bind_rows(choice_risky_gain, choice_risky_loss,choice_safe_gain,choice_safe_loss)
summary_choice <- choice_combined %>%
    group_by(Choice, Current) %>% 
  summarise(mean_resids = mean(mean_residuals,na.rm=TRUE),
            sem_resids = sd(mean_residuals, na.rm = TRUE) / sqrt(n()))

## `summarise()` has grouped output by 'Choice'. You can override using the
## `.groups` argument.

ggplot(summary_choice, aes(x = Choice, y = mean_resids,  fill=Current, group=Current)) +
    geom_bar(stat = "identity", position = position_dodge(0.7), width = 0.7, color = "black") +
    geom_errorbar(aes(ymin = mean_resids - sem_resids, ymax = mean_resids + sem_resids),
                  position = position_dodge(0.7), width = 0.2) + 
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
      scale_fill_manual(values = c("Gain" = "#FFD700", "Loss" = "#FF6347")) +

    labs(
      title = "Happiness Residuals (z-scored) Based on Choice",
      x = "Choice",
      y = "zHappy Resid"
    ) +
    theme_minimal() +
  ylim(-0.1, 0.1) + 
    poster_theme 

GAD SPLIT

#z-scored happiness
choice_risky_gain <- happyData1 %>%
  filter(Choice == 2 & Trial =="Gain" & !is.na(GAD_binary)) %>%
  group_by(userKey ,GAD_binary) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Current = "Gain")

choice_risky_loss <- happyData1 %>%
  filter(Choice == 2 & Trial =="Loss" & !is.na(GAD_binary)) %>%
  group_by(userKey,GAD_binary) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Current = "Loss")
  
# safe choice
choice_safe_gain <- happyData1 %>%
  filter(Choice == 1 & Trial=="Gain" & !is.na(GAD_binary)) %>%
  group_by(userKey,GAD_binary) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
        Current = "Gain")

choice_safe_loss <- happyData1 %>%
  filter(Choice == 1 & Trial=="Loss" & !is.na(GAD_binary)) %>%
  group_by(userKey,GAD_binary) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
        Current = "Loss")

  
choice_combined <- bind_rows(choice_risky_gain, choice_risky_loss,choice_safe_gain,choice_safe_loss)
summary_choice <- choice_combined %>%
    group_by(Choice,Current,GAD_binary) %>% 
  summarise(mean_happy = mean(mean_zhappy,na.rm=TRUE),
            mean_predhappy = mean(mean_zpredhappy,na.rm=TRUE),
            sem_happy = sd(mean_zhappy, na.rm = TRUE) / sqrt(n()),
            sem_predhappy = sd(mean_zpredhappy, na.rm = TRUE) / sqrt(n()))

## `summarise()` has grouped output by 'Choice', 'Current'. You can override using
## the `.groups` argument.

ggplot(summary_choice, aes(x = Choice, y = mean_happy,  color = GAD_binary)) +
    geom_point(position = position_dodge(0.5), width = 0.7,size=2) +
    geom_errorbar(aes(ymin = mean_happy - sem_happy, ymax = mean_happy + sem_happy),
                  position = position_dodge(0.5), width = 0.2) + 
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
  scale_color_manual(values = c("GAD 6+" = "#00BD9D", "GAD <6" = "#8BD7D2")) + 
    labs(
      title = "Happiness  (z-scored) Based on Choice",
      x = "Choice",
      y = "zHappy"
    ) + 
  facet_wrap(~Current) + 
    theme_minimal() +
    poster_theme 

## Warning in geom_point(position = position_dodge(0.5), width = 0.7, size = 2):
## Ignoring unknown parameters: `width`

Happiness Based on Choice and Future (Positive vs Negative)

happyData1 <- happyData1 %>%
  mutate(residuals=zHappy_filled - zHappyPred_filled)

#z-scored happiness
choice_risky_pos <- happyData1 %>%
  filter(Choice == 2 & NextIsland ==2 & LastIsland ==2 ) %>%
  group_by(userKey) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Future = "Positive")


choice_risky_neg <- happyData1 %>%
  filter(Choice == 2 & NextIsland ==1 & LastIsland ==1 ) %>%
  group_by(userKey) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Future = "Negative")
  
# safe choice
choice_safe_pos <- happyData1 %>%
  filter(Choice == 1 & NextIsland ==2 & LastIsland ==2 ) %>%
  group_by(userKey) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
         Future = "Positive")

choice_safe_neg <- happyData1 %>%
  filter(Choice == 1 & NextIsland ==1 & LastIsland ==1) %>%
  group_by(userKey) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
         Future = "Negative")

  
choice_combined <- bind_rows(choice_risky_pos, choice_risky_neg,choice_safe_pos,choice_safe_neg)
summary_choice <- choice_combined %>%
    group_by(Choice,Future) %>% 
  summarise(mean_happy = mean(mean_zhappy,na.rm=TRUE),
            mean_predhappy = mean(mean_zpredhappy,na.rm=TRUE),
            sem_happy = sd(mean_zhappy, na.rm = TRUE) / sqrt(n()),
            sem_predhappy = sd(mean_zpredhappy, na.rm = TRUE) / sqrt(n()))

## `summarise()` has grouped output by 'Choice'. You can override using the
## `.groups` argument.

ggplot(summary_choice, aes(x = Choice, y = mean_happy,  color = Choice)) +
    geom_point(position = position_dodge(1), width = 0.7) +
    geom_errorbar(aes(ymin = mean_happy - sem_happy, ymax = mean_happy + sem_happy),
                  position = position_dodge(1), width = 0.2, group="Choice") + 
      geom_point(aes(y = mean_predhappy), color = "lightblue", shape = 8, size = 3, position=position_dodge(0.7)) +
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    labs(
      title = "Happiness  (z-scored) Based on Choice",
      x = "Choice",
      y = "zHappy"
    ) +
   facet_wrap(~Future) +
    theme_minimal() +
    poster_theme 

## Warning in geom_point(position = position_dodge(1), width = 0.7): Ignoring
## unknown parameters: `width`

# residuals

choice_risky_pos <- happyData1 %>%
  filter(Choice == 2 & NextIsland ==2 & LastIsland ==2 ) %>%
  group_by(userKey) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Future = "Positive")

choice_risky_neg <- happyData1 %>%
  filter(Choice == 2 & NextIsland ==1 & LastIsland ==1 ) %>%
  group_by(userKey) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Future = "Negative")


  
# safe choice
choice_safe_pos <- happyData1 %>%
  filter(Choice == 1 & NextIsland ==2 & LastIsland ==2 ) %>%
  group_by(userKey) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
        Future = "Positive")

choice_safe_neg <- happyData1 %>%
  filter(Choice == 1  & NextIsland ==1 & LastIsland ==1) %>%
  group_by(userKey) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
        Future = "Negative")

  

choice_combined <- bind_rows(choice_risky_pos, choice_risky_neg,choice_safe_pos,choice_safe_neg)
summary_choice <- choice_combined %>%
    group_by(Choice, Future) %>% 
  summarise(mean_resids = mean(mean_residuals,na.rm=TRUE),
            sem_resids = sd(mean_residuals, na.rm = TRUE) / sqrt(n()))

## `summarise()` has grouped output by 'Choice'. You can override using the
## `.groups` argument.

ggplot(summary_choice, aes(x = Choice, y = mean_resids,  color = Choice)) +
    geom_point(position = position_dodge(1), width = 0.7) +
    geom_errorbar(aes(ymin = mean_resids - sem_resids, ymax = mean_resids + sem_resids),
                  position = position_dodge(1), width = 0.2, group="Choice") + 
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=
    labs(
      title = "Happiness Residuals (z-scored) Based on Choice",
      x = "Choice",
      y = "zHappy Resid"
    ) +
     facet_wrap(~Future) +
    theme_minimal() +
  ylim(-0.1, 0.1) + 
    poster_theme 

## Warning in geom_point(position = position_dodge(1), width = 0.7): Ignoring
## unknown parameters: `width`

Happiness Based on Choice, Future and Current Domain

happyData1 <- happyData1 %>%
  mutate(residuals = zHappy_filled - zHappyPred_filled)

# z-scored happiness with Trial filter
choice_risky_pos <- happyData1 %>%
  filter(Choice == 2 & NextIsland == 2 & LastIsland == 2 & (Trial == "Gain" | Trial == "Loss")) %>%
  group_by(userKey, Trial) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Future = "Positive")


choice_risky_neg <- happyData1 %>%
  filter(Choice == 2 & NextIsland == 1 & LastIsland == 1 & (Trial == "Gain" | Trial == "Loss")) %>%
  group_by(userKey, Trial) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Future = "Negative")

# Safe choice
choice_safe_pos <- happyData1 %>%
  filter(Choice == 1 & NextIsland == 2 & LastIsland == 2 & (Trial == "Gain" | Trial == "Loss")) %>%
  group_by(userKey, Trial) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
         Future = "Positive")

choice_safe_neg <- happyData1 %>%
  filter(Choice == 1 & NextIsland == 1 & LastIsland == 1 & (Trial == "Gain" | Trial == "Loss")) %>%
  group_by(userKey, Trial) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
         Future = "Negative")

# Combine all choices
choice_combined <- bind_rows(choice_risky_pos, choice_risky_neg, choice_safe_pos, choice_safe_neg)

# Summarise by Choice, Future, and Trial
summary_choice <- choice_combined %>%
  group_by(Choice, Future, Trial) %>%
  summarise(mean_happy = mean(mean_zhappy, na.rm = TRUE),
            mean_predhappy = mean(mean_zpredhappy, na.rm = TRUE),
            sem_happy = sd(mean_zhappy, na.rm = TRUE) / sqrt(n()),
            sem_predhappy = sd(mean_zpredhappy, na.rm = TRUE) / sqrt(n()))

## `summarise()` has grouped output by 'Choice', 'Future'. You can override using
## the `.groups` argument.

# Plot happiness based on choice, future, and trial
ggplot(summary_choice, aes(x = Choice, y = mean_happy, color = Trial)) +
  geom_point( width = 0.7) +
  geom_errorbar(aes(ymin = mean_happy - sem_happy, ymax = mean_happy + sem_happy),
                width = 0.2) + 
  geom_point(aes(y = mean_predhappy), color = "lightblue", shape = 8, size = 3) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
  labs(
    title = "Happiness (z-scored) Based on Choice and Trial",
    x = "Choice",
    y = "zHappy"
  ) +
  scale_color_manual(values = c("Gain" = "#FFD700", "Loss" = "#FF6347")) + 
  facet_wrap(~Future) +
  ylim(-0.2,0.2) +
  theme_minimal() +
  poster_theme

## Warning in geom_point(width = 0.7): Ignoring unknown parameters: `width`

# residuals

# Now for residuals with Trial filter
choice_risky_pos_res <- happyData1 %>%
  filter(Choice == 2 & NextIsland == 2 & LastIsland == 2 & (Trial == "Gain" | Trial == "Loss")) %>%
  group_by(userKey, Trial) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Future = "Positive")

choice_risky_neg_res <- happyData1 %>%
  filter(Choice == 2 & NextIsland == 1 & LastIsland == 1 & (Trial == "Gain" | Trial == "Loss")) %>%
  group_by(userKey, Trial) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Future = "Negative")


# Safe choice residuals
choice_safe_pos_res <- happyData1 %>%
  filter(Choice == 1 & NextIsland == 2 & LastIsland == 2 & (Trial == "Gain" | Trial == "Loss")) %>%
  group_by(userKey, Trial) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
         Future = "Positive")

choice_safe_neg_res <- happyData1 %>%
  filter(Choice == 1 & NextIsland == 1 & LastIsland == 1 & (Trial == "Gain" | Trial == "Loss")) %>%
  group_by(userKey, Trial) %>%
  summarise(mean_residuals = mean(residuals, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
         Future = "Negative")

# Combine all choices residuals
choice_combined_res <- bind_rows(choice_risky_pos_res, choice_risky_neg_res, choice_safe_pos_res, choice_safe_neg_res)

# Summarise by Choice, Future, and Trial for residuals
summary_choice_residuals <- choice_combined_res %>%
  group_by(Choice, Future, Trial) %>%
  summarise(mean_resids = mean(mean_residuals, na.rm = TRUE),
            sem_resids = sd(mean_residuals, na.rm = TRUE) / sqrt(n()))

## `summarise()` has grouped output by 'Choice', 'Future'. You can override using
## the `.groups` argument.

# Plot residuals based on choice, future, and trial
ggplot(summary_choice_residuals, aes(x = Choice, y = mean_resids, color = Trial)) +
  geom_point(position = position_dodge(1), width = 0.7) +
  geom_errorbar(aes(ymin = mean_resids - sem_resids, ymax = mean_resids + sem_resids),
                position = position_dodge(1), width = 0.2) +
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
  labs(
    title = "Happiness Residuals (z-scored) Based on Choice and Trial",
    x = "Choice",
    y = "zHappy Residuals"
  ) +
  facet_wrap(~Future) +
  theme_minimal() +
    scale_color_manual(values = c("Gain" = "#FFD700", "Loss" = "#FF6347")) + 
  ylim(-0.1, 0.1) + 
  poster_theme

## Warning in geom_point(position = position_dodge(1), width = 0.7): Ignoring
## unknown parameters: `width`

Split by GAD

happyData1 <- happyData1 %>%
  mutate(residuals = zHappy_filled - zHappyPred_filled)

# z-scored happiness with Trial filter
choice_risky_pos <- happyData1 %>%
  filter(Choice == 2 & NextIsland == 2 & LastIsland == 2 & (Trial == "Gain" | Trial == "Loss") & !is.na(GAD_binary)) %>%
  group_by(userKey, Trial, GAD_binary) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Future = "Positive")


choice_risky_neg <- happyData1 %>%
  filter(Choice == 2 & NextIsland == 1 & LastIsland == 1 & (Trial == "Gain" | Trial == "Loss") & !is.na(GAD_binary)) %>%
  group_by(userKey, Trial, GAD_binary) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Future = "Negative")

# Safe choice
choice_safe_pos <- happyData1 %>%
  filter(Choice == 1 & NextIsland == 2 & LastIsland == 2 & (Trial == "Gain" | Trial == "Loss") & !is.na(GAD_binary)) %>%
  group_by(userKey, Trial, GAD_binary) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
         Future = "Positive")

choice_safe_neg <- happyData1 %>%
  filter(Choice == 1 & NextIsland == 1 & LastIsland == 1 & (Trial == "Gain" | Trial == "Loss") & !is.na(GAD_binary)) %>%
  group_by(userKey, Trial, GAD_binary) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm = TRUE),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
         Future = "Negative")

# Combine all choices
choice_combined <- bind_rows(choice_risky_pos, choice_risky_neg, choice_safe_pos, choice_safe_neg)

# Summarise by Choice, Future, and Trial
summary_choice <- choice_combined %>%
  group_by(Choice, Future, Trial, GAD_binary) %>%
  summarise(mean_happy = mean(mean_zhappy, na.rm = TRUE),
            mean_predhappy = mean(mean_zpredhappy, na.rm = TRUE),
            sem_happy = sd(mean_zhappy, na.rm = TRUE) / sqrt(n()),
            sem_predhappy = sd(mean_zpredhappy, na.rm = TRUE) / sqrt(n()))

## `summarise()` has grouped output by 'Choice', 'Future', 'Trial'. You can
## override using the `.groups` argument.

# Plot happiness based on choice, future, and trial
ggplot(summary_choice, aes(x = Choice, y = mean_happy, color = GAD_binary)) +
  geom_point(position = position_dodge(0.5), size = 2) +
  geom_errorbar(aes(ymin = mean_happy - sem_happy, ymax = mean_happy + sem_happy),
                position = position_dodge(0.5), width = 0.2) + 
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
  labs(
    title = "Happiness (z-scored) Based on Choice, Trial, and GAD",
    x = "Choice",
    y = "zHappy"
  ) +
  scale_color_manual(values = c("GAD 6+" = "#00BD9D", "GAD <6" = "#8BD7D2")) + 
  facet_wrap(~ Future + Trial) +  # Facet by both Future and Trial
  ylim(-0.2, 0.2) +
  theme_minimal()

## GAD BINS

#z-scored happiness
choice_risky_gain <- happyData1 %>%
  filter(Choice == 2 & Trial =="Gain" & !is.na(GAD_bin)) %>%
  group_by(userKey ,GAD_bin) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
            GAD_score = first(gad7_total),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Current = "Gain")

cor.test(choice_risky_gain$GAD_score,choice_risky_gain$mean_zhappy,method="spearman")

## Warning in cor.test.default(choice_risky_gain$GAD_score,
## choice_risky_gain$mean_zhappy, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  choice_risky_gain$GAD_score and choice_risky_gain$mean_zhappy
## S = 1073789093, p-value = 0.04698
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.04569995

choice_risky_loss <- happyData1 %>%
  filter(Choice == 2 & Trial =="Loss" & !is.na(GAD_bin)) %>%
  group_by(userKey,GAD_bin) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
            GAD_score = first(gad7_total),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Current = "Loss")
  
cor.test(choice_risky_loss$GAD_score,choice_risky_loss$mean_zhappy,method="spearman")

## Warning in cor.test.default(choice_risky_loss$GAD_score,
## choice_risky_loss$mean_zhappy, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  choice_risky_loss$GAD_score and choice_risky_loss$mean_zhappy
## S = 1044594645, p-value = 0.1612
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.03282008

# safe choice
choice_safe_gain <- happyData1 %>%
  filter(Choice == 1 & Trial=="Gain" & !is.na(GAD_bin)) %>%
  group_by(userKey,GAD_bin) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
            GAD_score = first(gad7_total),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
        Current = "Gain")

cor.test(choice_safe_gain$GAD_score,choice_safe_gain$mean_zhappy,method="spearman")

## Warning in cor.test.default(choice_safe_gain$GAD_score,
## choice_safe_gain$mean_zhappy, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  choice_safe_gain$GAD_score and choice_safe_gain$mean_zhappy
## S = 1057276108, p-value = 0.4968
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.01576076

choice_safe_loss <- happyData1 %>%
  filter(Choice == 1 & Trial=="Loss" & !is.na(GAD_bin)) %>%
  group_by(userKey,GAD_bin) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
            GAD_score = first(gad7_total),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
        Current = "Loss")
cor.test(choice_safe_loss$GAD_score,choice_safe_loss$mean_zhappy,method="spearman")

## Warning in cor.test.default(choice_safe_loss$GAD_score,
## choice_safe_loss$mean_zhappy, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  choice_safe_loss$GAD_score and choice_safe_loss$mean_zhappy
## S = 1104933419, p-value = 0.5979
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.01219951

choice_combined <- bind_rows(choice_risky_gain, choice_risky_loss,choice_safe_gain,choice_safe_loss)
summary_choice <- choice_combined %>%
    group_by(Choice,Current,GAD_bin) %>% 
  summarise(mean_happy = mean(mean_zhappy,na.rm=TRUE),
            mean_predhappy = mean(mean_zpredhappy,na.rm=TRUE),
            sem_happy = sd(mean_zhappy, na.rm = TRUE) / sqrt(n()),
            sem_predhappy = sd(mean_zpredhappy, na.rm = TRUE) / sqrt(n()))

## `summarise()` has grouped output by 'Choice', 'Current'. You can override using
## the `.groups` argument.

ggplot(summary_choice, aes(x = Choice, y = mean_happy,  color = GAD_bin)) +
    geom_point(position = position_dodge(0.5), width = 0.7,size=2) +
    geom_errorbar(aes(ymin = mean_happy - sem_happy, ymax = mean_happy + sem_happy),
                  position = position_dodge(0.5), width = 0.2) + 
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
    scale_color_manual(values = c("0-4" = "#8BD7D2", "5-9" = "#6FACA8", "10-14" = "#53817E", "15-21" = "#385654")) +
    labs(
      title = "Happiness  (z-scored) Based on Choice",
      x = "Choice",
      y = "zHappy"
    ) + 
  facet_wrap(~Current) + 
    theme_minimal() +
    poster_theme

## Warning in geom_point(position = position_dodge(0.5), width = 0.7, size = 2):
## Ignoring unknown parameters: `width`

## PHQ bins

#z-scored happiness
choice_risky_gain <- happyData1 %>%
  filter(Choice == 2 & Trial =="Gain" & !is.na(PHQ_bin)) %>%
  group_by(userKey ,PHQ_bin) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
            PHQ_score = first(phq8_total),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Current = "Gain")

cor.test(choice_risky_gain$PHQ_score,choice_risky_gain$mean_zhappy,method="spearman")

## Warning in cor.test.default(choice_risky_gain$PHQ_score,
## choice_risky_gain$mean_zhappy, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  choice_risky_gain$PHQ_score and choice_risky_gain$mean_zhappy
## S = 1129183255, p-value = 0.314
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.02306938

choice_risky_loss <- happyData1 %>%
  filter(Choice == 2 & Trial =="Loss" & !is.na(PHQ_bin)) %>%
  group_by(userKey,PHQ_bin) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
            PHQ_score = first(phq8_total),
              .groups = 'drop') %>%
  mutate(Choice = "Risky",
         Current = "Loss")
  
cor.test(choice_risky_loss$PHQ_score,choice_risky_loss$mean_zhappy,method="spearman")

## Warning in cor.test.default(choice_risky_loss$PHQ_score,
## choice_risky_loss$mean_zhappy, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  choice_risky_loss$PHQ_score and choice_risky_loss$mean_zhappy
## S = 1074828825, p-value = 0.1309
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.03523087

# safe choice
choice_safe_gain <- happyData1 %>%
  filter(Choice == 1 & Trial=="Gain" & !is.na(PHQ_bin)) %>%
  group_by(userKey,PHQ_bin) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
            PHQ_score = first(phq8_total),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
        Current = "Gain")

cor.test(choice_safe_gain$PHQ_score,choice_safe_gain$mean_zhappy,method="spearman")

## Warning in cor.test.default(choice_safe_gain$PHQ_score,
## choice_safe_gain$mean_zhappy, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  choice_safe_gain$PHQ_score and choice_safe_gain$mean_zhappy
## S = 1060923064, p-value = 0.09156
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.03894443

choice_safe_loss <- happyData1 %>%
  filter(Choice == 1 & Trial=="Loss" & !is.na(PHQ_bin)) %>%
  group_by(userKey,PHQ_bin) %>%
  summarise(mean_zhappy = mean(zHappy_filled, na.rm = TRUE),
            mean_zpredhappy = mean(zHappyPred_filled, na.rm=TRUE),
            PHQ_score = first(phq8_total),
              .groups = 'drop') %>%
  mutate(Choice = "Safe",
        Current = "Loss")
cor.test(choice_safe_loss$PHQ_score,choice_safe_loss$mean_zhappy,method="spearman")

## Warning in cor.test.default(choice_safe_loss$PHQ_score,
## choice_safe_loss$mean_zhappy, : Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  choice_safe_loss$PHQ_score and choice_safe_loss$mean_zhappy
## S = 1123559744, p-value = 0.9408
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##          rho 
## -0.001708993

choice_combined <- bind_rows(choice_risky_gain, choice_risky_loss,choice_safe_gain,choice_safe_loss)
summary_choice <- choice_combined %>%
    group_by(Choice,Current,PHQ_bin) %>% 
  summarise(mean_happy = mean(mean_zhappy,na.rm=TRUE),
            mean_predhappy = mean(mean_zpredhappy,na.rm=TRUE),
            sem_happy = sd(mean_zhappy, na.rm = TRUE) / sqrt(n()),
            sem_predhappy = sd(mean_zpredhappy, na.rm = TRUE) / sqrt(n()))

## `summarise()` has grouped output by 'Choice', 'Current'. You can override using
## the `.groups` argument.

ggplot(summary_choice, aes(x = Choice, y = mean_happy,  color = PHQ_bin)) +
    geom_point(position = position_dodge(0.5), width = 0.7,size=2) +
    geom_errorbar(aes(ymin = mean_happy - sem_happy, ymax = mean_happy + sem_happy),
                  position = position_dodge(0.5), width = 0.2) + 
    geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
      scale_color_manual(values = c("0-4" = "#FFC0CB", "5-9" = "#CC9AA2", "10-14" = "#99737A", "15-24" = "#664D51")) +
    labs(
      title = "Happiness  (z-scored) Based on Choice",
      x = "Choice",
      y = "zHappy"
    ) + 
  facet_wrap(~Current) + 
    theme_minimal() +
    poster_theme

## Warning in geom_point(position = position_dodge(0.5), width = 0.7, size = 2):
## Ignoring unknown parameters: `width`

happyData1_ev <- happyData1 %>%
  mutate(SafeEV = SafeProb*SafeValue,
         RiskyEV = RiskyProb*RiskyValue)

risky_data <- happyData1_ev %>%
  group_by(userKey, Trial, RiskyEV) %>%
  summarise(
    total_choices = n(),
    risky_choices = sum(Choice == 2, na.rm = TRUE),  # Assuming Choice == 2 indicates a risky choice
    percent_risky = (risky_choices / total_choices) * 100,
    .groups = 'drop'
  )

# Calculate the mean percent_risky per participant for each abs(RiskyEV)
mean_risky_ev_data <- risky_data %>%
  mutate(abs_RiskyEV = abs(RiskyEV)) %>%         # Use absolute RiskyEV
  group_by(userKey, Trial, abs_RiskyEV) %>%
  summarise(mean_percent_risky = mean(percent_risky, na.rm = TRUE),
            .groups = 'drop')

# Now calculate the average across participants for each abs(RiskyEV), split by Trial
overall_mean_risky_ev_data <- mean_risky_ev_data %>%
  group_by(Trial, abs_RiskyEV) %>%
  summarise(mean_percent_risky_all = mean(mean_percent_risky, na.rm = TRUE),
            .groups = 'drop')

# Plot the data
ggplot(overall_mean_risky_ev_data, aes(x = abs_RiskyEV, y = mean_percent_risky_all, color = Trial)) +
  geom_line(size = 1) +
  geom_hline(yintercept = 50, linetype = "dotted", color = "black") +  # Add a dotted line at y=0
  geom_point() +                                 # Add points for clarity
  labs(x = "Absolute RiskyEV", y = "Mean Percentage of Risky Choices",
       title = "Mean % of Risky Choices by Absolute RiskyEV",
       color = "Trial") +
  theme_minimal() +
  ylim(0, 100) +
  scale_color_manual(values = c("Gain" = "#FFD700", "Loss" = "#FF6347"))  +
  poster_theme

## split by GAD
risky_data <- happyData1_ev %>%
  filter(!is.na(GAD_binary)) %>%
  group_by(userKey, Trial, RiskyEV, GAD_binary) %>%
  summarise(
    total_choices = n(),
    risky_choices = sum(Choice == 2, na.rm = TRUE),  # Assuming Choice == 2 indicates a risky choice
    percent_risky = (risky_choices / total_choices) * 100,
    .groups = 'drop'
  )

mean_risky_ev_data_gad <- risky_data %>%
  mutate(abs_RiskyEV = abs(RiskyEV)) %>%
  group_by(userKey, Trial, GAD_binary, abs_RiskyEV) %>%
  summarise(mean_percent_risky = mean(percent_risky, na.rm = TRUE),
            .groups = 'drop')

overall_mean_risky_ev_data_gad <- mean_risky_ev_data_gad %>%
  group_by(Trial, GAD_binary, abs_RiskyEV) %>%
  summarise(mean_percent_risky_all = mean(mean_percent_risky, na.rm = TRUE),
            .groups = 'drop')

# Create the plot
ggplot(overall_mean_risky_ev_data_gad, aes(x = abs_RiskyEV, y = mean_percent_risky_all, color = GAD_binary, linetype = Trial)) +
  geom_line(size = 1) +   
    geom_hline(yintercept = 50, linetype = "dotted", color = "black") + 
  geom_point() +                                              # Points for clarity
  labs(x = "Absolute RiskyEV", y = "Mean Percentage of Risky Choices",
       title = "Mean % of Risky Choices by Absolute RiskyEV (split by GAD)",
       color = "GAD Group", linetype = "Trial Type") +
  theme_minimal() +
  ylim(0, 100) +
  scale_color_manual(values = c("GAD 6+" = "#00BD9D", "GAD <6" = "#8BD7D2")) +  # Custom colors for GAD groups
  scale_linetype_manual(values = c("Gain" = "solid", "Loss" = "dotted")) + 
  poster_theme 

  ## PHQ

risky_data <- happyData1_ev %>%
  filter(!is.na(PHQ_binary)) %>%
  group_by(userKey, Trial, RiskyEV, PHQ_binary) %>%
  summarise(
    total_choices = n(),
    risky_choices = sum(Choice == 2, na.rm = TRUE),  # Assuming Choice == 2 indicates a risky choice
    percent_risky = (risky_choices / total_choices) * 100,
    .groups = 'drop'
  )

mean_risky_ev_data_phq <- risky_data %>%
  mutate(abs_RiskyEV = abs(RiskyEV)) %>%
  group_by(userKey, Trial, PHQ_binary, abs_RiskyEV) %>%
  summarise(mean_percent_risky = mean(percent_risky, na.rm = TRUE),
            .groups = 'drop')

overall_mean_risky_ev_data_phq <- mean_risky_ev_data_phq %>%
  group_by(Trial, PHQ_binary, abs_RiskyEV) %>%
  summarise(mean_percent_risky_all = mean(mean_percent_risky, na.rm = TRUE),
            .groups = 'drop')

# Create the plot
ggplot(overall_mean_risky_ev_data_phq, aes(x = abs_RiskyEV, y = mean_percent_risky_all, color = PHQ_binary, linetype = Trial)) +
  geom_line(size = 1) +   
    geom_hline(yintercept = 50, linetype = "dotted", color = "black") + 
  geom_point() +                                              # Points for clarity
  labs(x = "Absolute RiskyEV", y = "Mean Percentage of Risky Choices",
       title = "Mean % of Risky Choices by Absolute RiskyEV (split by PHQ)",
       color = "PHQ Group", linetype = "Trial Type") +
  theme_minimal() +
  ylim(0, 100) +
  scale_color_manual(values = c("PHQ 7+" = "red", "PHQ <7" = "pink")) +  # Custom colors for GAD groups
  scale_linetype_manual(values = c("Gain" = "solid", "Loss" = "dotted")) + 
  poster_theme 

happyData1_ev <- happyData1_ev %>%
  mutate(RiskyEV_bin = cut(abs(RiskyEV), 
                           breaks = c(0, 9, 19, 29, 39, 49, Inf), 
                           labels = c("0-9", "10-19", "20-29", "30-39", "40-49", "50+"),
                           right = TRUE))

happy_data <- happyData1_ev %>%
  filter(!is.na(GAD_binary)) %>%
  group_by(userKey, Trial, RiskyEV_bin,GAD_binary) %>%
  summarise(
    mean_zhappy = mean(zHappy,na.rm=TRUE)
  ) 

## `summarise()` has grouped output by 'userKey', 'Trial', 'RiskyEV_bin'. You can
## override using the `.groups` argument.

# Calculate the mean percent_risky per participant for each abs(RiskyEV)
mean_happy_ev_data <- happy_data %>%
  mutate(abs_RiskyEV = RiskyEV_bin) %>%         # Use absolute RiskyEV
  group_by(userKey, Trial, abs_RiskyEV) %>%
  summarise(mean_z_happy = mean(mean_zhappy, na.rm = TRUE),
            .groups = 'drop')

# Now calculate the average across participants for each abs(RiskyEV), split by Trial
overall_mean_happy_ev_data <- mean_happy_ev_data %>%
  group_by(Trial, abs_RiskyEV) %>%
  summarise(mean_z_happy_all = mean(mean_z_happy, na.rm = TRUE),
            sem_z_happy_all = sd(mean_z_happy, na.rm=TRUE) / sqrt(n()),
            .groups = 'drop')

# Plot the data with error bars
ggplot(overall_mean_happy_ev_data, aes(x = abs_RiskyEV, y = mean_z_happy_all, color = Trial, group = Trial)) +
  geom_line(size = 1) +                             # Ensure the line connects the points
  geom_point() +                                    # Points for clarity
  geom_errorbar(aes(ymin = mean_z_happy_all - sem_z_happy_all, ymax = mean_z_happy_all + sem_z_happy_all), width = 0.2) +  # Error bars for SEM
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Dotted line at y=0
  labs(x = "Absolute RiskyEV", y = "Mean z Happy Score",
       title = "Mean z Happy Score by Absolute RiskyEV",
       color = "Trial") +
  theme_minimal() +
  scale_color_manual(values = c("Gain" = "#FFD700", "Loss" = "#FF6347")) +  # Custom colors
  poster_theme

## GAD Split


# Calculate the mean z_happy per participant for each abs(RiskyEV), including GAD groupings
mean_happy_ev_data_gad <- happy_data %>%
  mutate(abs_RiskyEV = RiskyEV_bin) %>%         # Use absolute RiskyEV
  group_by(userKey, Trial, GAD_binary, abs_RiskyEV) %>%
  summarise(mean_z_happy = mean(mean_zhappy, na.rm = TRUE),
            .groups = 'drop')

# Now calculate the average across participants for each abs(RiskyEV), split by Trial and GAD group
overall_mean_happy_ev_data_gad <- mean_happy_ev_data_gad %>%
  group_by(Trial, GAD_binary, abs_RiskyEV) %>%
  summarise(mean_z_happy_all = mean(mean_z_happy, na.rm = TRUE),
             sem_z_happy_all = sd(mean_z_happy, na.rm=TRUE) / sqrt(n()),
            .groups = 'drop')

# Plot the data with GAD splitting
ggplot(overall_mean_happy_ev_data_gad, aes(x = abs_RiskyEV, y = mean_z_happy_all, color = GAD_binary, linetype = Trial, group = interaction(GAD_binary, Trial))) +
  geom_line(size = 1) +                             # Ensure the line connects the points
  geom_point() +                                    # Points for clarity
  geom_errorbar(aes(ymin = mean_z_happy_all - sem_z_happy_all, ymax = mean_z_happy_all + sem_z_happy_all), width = 0.2) +  # Error bars for SEM
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Dotted line at y=0
  labs(x = "Absolute RiskyEV", y = "Mean z Happy Score",
       title = "Mean z Happy Score by Absolute RiskyEV (split by GAD)",
       color = "GAD Group", linetype = "Trial Type") +
  theme_minimal() +
  scale_color_manual(values = c("GAD 6+" = "#00BD9D", "GAD <6" = "#8BD7D2")) +  # Custom colors for GAD groups
  scale_linetype_manual(values = c("Gain" = "solid", "Loss" = "dotted")) +      # Solid for Gain, dotted for Loss
  poster_theme

# Calculate the mean z_happy per participant for each abs(RiskyEV), including GAD groupings
happy_data <- happyData1_ev %>%
  filter(!is.na(PHQ_binary)) %>%
  group_by(userKey, Trial, RiskyEV_bin,PHQ_binary) %>%
  summarise(
    mean_zhappy = mean(zHappy,na.rm=TRUE)
  ) 

## `summarise()` has grouped output by 'userKey', 'Trial', 'RiskyEV_bin'. You can
## override using the `.groups` argument.

mean_happy_ev_data_phq <- happy_data %>%
  mutate(abs_RiskyEV = RiskyEV_bin) %>%         # Use absolute RiskyEV
  group_by(userKey, Trial, PHQ_binary, abs_RiskyEV) %>%
  summarise(mean_z_happy = mean(mean_zhappy, na.rm = TRUE),
            .groups = 'drop')

# Now calculate the average across participants for each abs(RiskyEV), split by Trial and GAD group
overall_mean_happy_ev_data_phq <- mean_happy_ev_data_phq %>%
  group_by(Trial, PHQ_binary, abs_RiskyEV) %>%
  summarise(mean_z_happy_all = mean(mean_z_happy, na.rm = TRUE),
             sem_z_happy_all = sd(mean_z_happy, na.rm=TRUE) / sqrt(n()),
            .groups = 'drop')

# Plot the data with GAD splitting
ggplot(overall_mean_happy_ev_data_phq, aes(x = abs_RiskyEV, y = mean_z_happy_all, color = PHQ_binary, linetype = Trial, group = interaction(PHQ_binary, Trial))) +
  geom_line(size = 1) +                             # Ensure the line connects the points
  geom_point() +                                    # Points for clarity
  geom_errorbar(aes(ymin = mean_z_happy_all - sem_z_happy_all, ymax = mean_z_happy_all + sem_z_happy_all), width = 0.2) +  # Error bars for SEM
  geom_hline(yintercept = 0, linetype = "dotted", color = "black") +  # Dotted line at y=0
  labs(x = "Absolute RiskyEV", y = "Mean z Happy Score",
       title = "Mean z Happy Score by Absolute RiskyEV (split by PHQ)",
       color = "PHQ Group", linetype = "Trial Type") +
  theme_minimal() +
  scale_color_manual(values = c("PHQ 7+" = "red", "PHQ <7" = "pink")) +  # Custom colors for GAD groups
  scale_linetype_manual(values = c("Gain" = "solid", "Loss" = "dotted")) +      # Solid for Gain, dotted for Loss
  poster_theme

Gender x Anxiety Analyses

T_initquiz <- process_initquiz(initquiz_file)
happyData1_demo <- merge(happyData1, T_initquiz, by="userKey")

length(unique(happyData1_demo$userKey)) #1429

## [1] 1429

anxiety, gender and z-scored happiness

risky_data <- happyData1_demo %>%
    group_by(userKey, Trial,gender_group) %>%
    summarise(total_choices = n(),
              risky_choices = sum(Choice == 2, na.rm = TRUE),
              percent_risky = (risky_choices / total_choices) * 100,
              GAD_total = mean(gad7_total, na.rm = TRUE),  # Add GAD scores
             PHQ_total = mean(phq8_total, na.rm = TRUE),
              .groups = 'drop')


summary_data <- risky_data %>%
    group_by(Trial,gender_group) %>%
    summarise(mean_percent_risky = mean(percent_risky, na.rm = TRUE),
              se_percent_risky = sd(percent_risky, na.rm = TRUE) / sqrt(n()),
              .groups = 'drop') %>%
  na.omit()
  
  ggplot(summary_data, aes(x = Trial, y = mean_percent_risky, color = gender_group)) +
    geom_point(position = position_dodge(0.8), size = 3) +
    geom_errorbar(aes(ymin = mean_percent_risky - se_percent_risky, ymax = mean_percent_risky + se_percent_risky),position=position_dodge(0.8),
                  width = 0.2) +
    labs(
      title = "% Risky Choices (Play 1) N=1429",
      x = "Trial Type",
      y = "% Risky Choices"
    ) +
    theme_minimal() +
    ylim(40, 60) +  
    poster_theme + 
    scale_color_manual(values = c("Male" = "lightblue", "Female" = "pink","Other" = "purple"))

# correlations with GAD and PHQ
risky_data_males <- risky_data %>%
  filter(gender_group == "Male") %>%
  na.omit()

risky_data_females <- risky_data %>%
  filter(gender_group == "Female") %>%
  na.omit()

# males gain and loss gad
cor.test(risky_data_males$percent_risky[risky_data_males$Trial == "Gain"], risky_data_males$GAD_total[risky_data_males$Trial == "Gain"],method="spearman")

## Warning in
## cor.test.default(risky_data_males$percent_risky[risky_data_males$Trial == :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  risky_data_males$percent_risky[risky_data_males$Trial == "Gain"] and risky_data_males$GAD_total[risky_data_males$Trial == "Gain"]
## S = 3743136, p-value = 0.879
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## 0.009092746

 cor.test(risky_data_males$percent_risky[risky_data_males$Trial == "Loss"], risky_data_males$GAD_total[risky_data_males$Trial == "Loss"],method="spearman")

## Warning in
## cor.test.default(risky_data_males$percent_risky[risky_data_males$Trial == :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  risky_data_males$percent_risky[risky_data_males$Trial == "Loss"] and risky_data_males$GAD_total[risky_data_males$Trial == "Loss"]
## S = 3594234, p-value = 0.4162
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.04851107

# males gain and loss phq
cor.test(risky_data_males$percent_risky[risky_data_males$Trial == "Gain"], risky_data_males$PHQ_total[risky_data_males$Trial == "Gain"],method="spearman")

## Warning in
## cor.test.default(risky_data_males$percent_risky[risky_data_males$Trial == :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  risky_data_males$percent_risky[risky_data_males$Trial == "Gain"] and risky_data_males$PHQ_total[risky_data_males$Trial == "Gain"]
## S = 3757944, p-value = 0.931
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## 0.005172813

cor.test(risky_data_males$percent_risky[risky_data_males$Trial == "Loss"], risky_data_males$PHQ_total[risky_data_males$Trial == "Loss"],method="spearman")

## Warning in
## cor.test.default(risky_data_males$percent_risky[risky_data_males$Trial == :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  risky_data_males$percent_risky[risky_data_males$Trial == "Loss"] and risky_data_males$PHQ_total[risky_data_males$Trial == "Loss"]
## S = 3483311, p-value = 0.1915
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.07787535

# females gain and loss gad
cor.test(risky_data_females$percent_risky[risky_data_females$Trial == "Gain"], risky_data_females$GAD_total[risky_data_females$Trial == "Gain"],method="spearman")

## Warning in
## cor.test.default(risky_data_females$percent_risky[risky_data_females$Trial == :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  risky_data_females$percent_risky[risky_data_females$Trial == "Gain"] and risky_data_females$GAD_total[risky_data_females$Trial == "Gain"]
## S = 196509047, p-value = 0.2198
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.03753965

cor.test(risky_data_females$percent_risky[risky_data_females$Trial == "Loss"], risky_data_females$GAD_total[risky_data_females$Trial == "Loss"],method="spearman")

## Warning in
## cor.test.default(risky_data_females$percent_risky[risky_data_females$Trial == :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  risky_data_females$percent_risky[risky_data_females$Trial == "Loss"] and risky_data_females$GAD_total[risky_data_females$Trial == "Loss"]
## S = 207412250, p-value = 0.6043
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##         rho 
## -0.01586196

# females gain and loss phq
cor.test(risky_data_females$percent_risky[risky_data_females$Trial == "Gain"], risky_data_females$PHQ_total[risky_data_females$Trial == "Gain"],method="spearman")

## Warning in
## cor.test.default(risky_data_females$percent_risky[risky_data_females$Trial == :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  risky_data_females$percent_risky[risky_data_females$Trial == "Gain"] and risky_data_females$PHQ_total[risky_data_females$Trial == "Gain"]
## S = 205682633, p-value = 0.8092
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##          rho 
## -0.007390658

cor.test(risky_data_females$percent_risky[risky_data_females$Trial == "Loss"], risky_data_females$PHQ_total[risky_data_females$Trial == "Loss"],method="spearman") # no longer correlation between risk taking and PHQ

## Warning in
## cor.test.default(risky_data_females$percent_risky[risky_data_females$Trial == :
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  risky_data_females$percent_risky[risky_data_females$Trial == "Loss"] and risky_data_females$PHQ_total[risky_data_females$Trial == "Loss"]
## S = 197861251, p-value = 0.3123
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##        rho 
## 0.03091684